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Summary 


The theory of operation of spring gravity meters on a stabilized 
platform in the presence of horizontal and vertical accelerations is con- 
sidered in the first section of this paper. It is shown that levelling 
errors with the same period as the horizontal accelerations are very 
important. ‘The amplitude of such levelling errors should not exceed 
about 8 seconds of arc when the amplitude of the horizontal accelera- 
tions is 50000 mgal. A second effect, which is present even on a per- 
fectly stabilized platform, has been named the ‘‘cross-coupling”’ effect. 
It is produced by coupling between the vertical and horizontal 
accelerations and may, for a typical gravity meter, cause an error in 
excess of 100 mgal when the horizontal and vertical accelerations both 
have amplitudes of 50000 mgal. 

The second part of the paper is devoted to a discussion of a gravity 
meter free to swing in gimbals. An expression for the second-order 
correction is derived taking into account free and forced oscillations 
of the gimbal system. This expression takes an especially simple 
form 482 (0 is the deflection of the gimbal system from the true 
vertical) when the gravity sensing element is positioned at the correct 
distance below the gimbal axis. 


1. Introduction 


Several surface ship gravity meters have been developed in the last few years 
and more are in process of development. It is therefore important that the theory 


of operation of gravity meters in the presence of large horizontal and vertical 
accelerations should be understood. 
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Summary 

The theory of operation of spring gravity meters on a stabilized 
platform in the presence of horizontal and vertical accelerations is con- 
sidered in the first section of this paper. It is shown that levelling 
errors with the same period as the horizontal accelerations are very 
important. The amplitude of such levelling errors should not exceed 
about 8 seconds of arc when the amplitude of the horizontal accelera- 
tions is 50000 mgal. A second effect, which is present even on a per- 
fectly stabilized platform, has been named the “‘cross-coupling”’ effect. 
It is produced by coupling between the vertical and horizontal 
accelerations and may, for a typical gravity meter, cause an error in 
excess of 100 mgal when the horizontal and vertical accelerations both 
have amplitudes of 50000 mgal. 

The second part of the paper is devoted to a discussion of a gravity 
meter free to swing in gimbals. An expression for the second-order 
correction is derived taking into account free and forced oscillations 
of the gimbal system. This expression takes an especially simple 
form }g02 (0 is the deflection of the gimbal system from the true 


vertical) when the gravity sensing element is positioned at the correct 
distance below the gimbal axis. 


1. Introduction 


Several surface ship gravity meters have been developed in the last few years 
and more are in process of development. It is therefore important that the theory 
of operation of gravity meters in the presence of large horizontal and vertical 
accelerations should be understood. 
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Two distinct methods may be used to eliminate the effect of horizontal 
accelerations. The meter may be stabilized to hold its measuring axis aligned 
along the true vertical so that only the vertical component of acceleration is 
measured, or the meter may be allowed to swing freely in gimbals. In the latter 
case the total resultant acceleration, which is always larger than the vertical 
component, is measured and a correction, the second-order or Browne correction, 
must be subtracted. 


2. Theory of gravity meters on a stabilized platform 

A gravity meter for use on a stabilized platform must be constructed to be 
unaffected by accelerations acting at right angles to its measuring axis. Consider 
such a meter on an almost horizontal platform such that its measuring axis makes 
angles « and 8 with two perpendicular horizontal axes X and Y and an angle y 


Stabilized 
platform 


Measuring axis of 
gravity meter 


Vertical 


Fic. 1.—Gravity meter on a stable platform. 


with the vertical Z (see Figure 1). Put a = 7/2—« and b = /2—8, so that a, b 
and y are small, and 
sin*a + sin*b + cos*y = 1. (1) 


Third- and higher-order terms in a, b, y will be neglected. In the absence of disturb- 
ing accelerations the acceleration measured, A, is 


A = goosy = (2) 
If g—A = 10~%g, i.e. the measuring error is 1 mgal, 
= 2. radians? 


y=5'. 
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On a stationary platform a levelling precision of only 5’ is needed to obtain 
1 mgal accuracy. 


Now let the meter and table be subjected to periodic accelerations #, j and 2 
where 


# = sin wt 
= Hrsin(wet + 8) (3) 
= sin(wst +) 
and consider a and 6 to vary as 
@ = +A) 
b = bo +b; sin(vet + 4). 
The measured acceleration, A, is given by 
= —Xcosa—ycosB+(g+2)cosy 
= sin wit sin{ao + a; sin(v;t +A)} 
— Sin(wot + 8) sin{bo + 5; sin(vet + ¢)} 
+ €)} x —sin®{ap + a) sin(vit + A)} — 
—sin?{bo + sin(vet + $}}]. (5) 


It is beyond the scope of this paper to try to treat the preceding equation 
rigorously. However since we are interested in the case in which ao, bo, a; and 5; 
are small compared to 1, second- and higher-order terms in the expansion of 
equation (5) will be neglected. Making this expansion of equation (5) and taking 
a time average over a time long compared with any of the periods concerned, 


A-g = Jar? + 
— a, sin wt + A) 
— by; sin(wot + 8) sin(vet + 
— sin(wst + €){aga; sin(vyt + A) + bob; sin(vet + (6) 
Of these terms, the first is the equivalent of that occurring in the stationary case 
and the others have a zero time average unless the periods of the sine terms are 
equal. With w; = »;, the time average of the second term is — }a)#, cos A and 
with w2 = ve, the time average of the third is — $5; cos(—¢). These two terms 
depend on the first power of the levelling error and are therefore more important 


than the first and fourth which are of the second order in this small quantity. 
The error introduced by the second term is 


cos A. (7) 


In the worst case where A = 0, 7, ... etc., so that cosA = +1 the magnitude of 
the error is 


$a 


If |A—g| < 1mgal and #; = 50000mgal then a; must be less than 8”. Thus 
stabilization errors with the same period as the horizontal accelerations are 
extremely important. 
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Two distinct methods may be used to eliminate the effect of horizontal 
accelerations. The meter may be stabilized to hold its measuring axis aligned 
along the true vertical so that only the vertical component of acceleration is 
measured, or the meter may be allowed to swing freely in gimbals. In the latter 
case the total resultant acceleration, which is always larger than the vertical 
component, is measured and a correction, the second-order or Browne correction, 
must be subtracted. 


2. Theory of gravity meters on a stabilized platform 

A gravity meter for use on a stabilized platform must be constructed to be 
unaffected by accelerations acting at right angles to its measuring axis. Consider 
such a meter on an almost horizontal platform such that its measuring axis makes 
angles « and 8 with two perpendicular horizontal axes X and Y and an angle y 


stabilized 
platform 


Measuring axis of 
gravity meter 


Vertical 


Fic. 1.—Gravity meter on a stable platform. 


with the vertical Z (see Figure 1). Put a = 7/2— and b = 2/2—8, so that a, b 
and y are small, and 


sin?a + sin*b + cos*y = 1. (1) 


Third- and higher-order terms in a, b, y will be neglected. In the absence of disturb- 
ing accelerations the acceleration measured, A, is 


A = geosy = g(1—}/*). (2) 
If g—A = 10~g, i.e. the measuring error is 1 mgal, 
y? = 2. radians? 


y= 
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On a stationary platform a levelling precision of only 5’ is needed to obtain 
1 mgal accuracy. 


Now let the meter and table be subjected to periodic accelerations #, j and # 
where 
= Xsinwt 
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It is beyond the scope of this paper to try to treat the preceding equation 
rigorously. However since we are interested in the case in which ao, bo, a; and 6; 
are small compared to 1, second- and higher-order terms in the expansion of 
equation (5) will be neglected. Making this expansion of equation (5) and taking 
a time average over a time long compared with any of the periods concerned, 


Ang = + + far? + 
— sin sin(v;t + A) 
— sin(wet + 8) sin(vet +4) 
— sin(wgt + €){aoa) sin(vyt + A) + bob; sin(vet + ¢4)}. (6) 


Of these terms, the first is the equivalent of that occurring in the stationary case 
and the others have a zero time average unless the periods of the sine terms are 
equal. With w; = », the time average of the second term is — }a)#; cos A and 
with we = ve, the time average of the third is — 4417; cos(8—¢). These two terms 
depend on the first power of the levelling error and are therefore more important 
than the first and fourth which are of the second order in this small quantity. 
The error introduced by the second term is 


cos A. (7) 


In the worst case where A = 0, z,... etc., so that cosA = +1 the magnitude of 
the error is 


If \A —g| < 1mgal and #; = 50000mgal then a; must be less than 8”. Thus 
stabilization errors with the same period as the horizontal accelerations are 
extremely important. 
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Suppose the gravity meter to be mounted in a gimbal system with a natural 
period long with respect to the wave periods, or to be mcunted on a stabilized 
platform controlled by a long-period pendulum. Suppose jj = o and let the 
undamped free period of the gimbal system or controlling pendulum be 
To = 2m/wo? where wp < wm. The periodic deflection, a, of the gravity meter’s 
axis produced by the acceleration #,, neglecting friction in the gimbal system or 
controlling pendulum, is given by the equation 


whose solution for wo < w is 


wort, 
sin (8) 
wg 


Comparing the expressions for a in (4) and (8), 


a, = Wo?/w;? . ¥1/g, vy = wanddA =o. 
By (7), the error introduced into the gravity measurement is 
. (9) 


The minimum natural period that will reduce this error below 1 mgal for a 
given wave frequency w can be estimated from the inequality 


2 
Wo? < 2. 1065 


Consider the example #/g = 0-05 as before with a 10 s wave period 
wo < 2. 
To > 103/24/2 = 3548. 


Hence the free period of the gimbal suspension or controlling reference pendulum 
must be at least 6 minutes to obtain 1 mgal accuracy under the above conditions. 

Two types of stable elements are used for controlling stabilized platforms. 
A Schuler-tuned inertial system with its free period of 84 minutes can probably 
meet the demands required of a platform on which to make gravity measurements 
at sea. A more common stable element consists of a gyro (or gyros) controlled 
by a pendulum or liquid level, so that the gyro is made to precess towards the 
vertical indicated by the controlling device. Such an element will suffer from 
periodic errors of the same periods as the horizontal accelerations. 

A second source of error in the stabilized platform is introduced by the servo 
motors which cause the platform to remain aligned with the stable element. Servo 
systems suffer from errors proportional to the velocity and acceleration of the 
motion which they follow and to the load impressed on them. All three of these 
quantities contain components with periods equal to those of the motions of the 
ship. Very accurate servos can be made but they are hard to keep operating 
perfectly, and it is difficult to tell when they are not operating quite correctly. 
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The cross-coupling effect is caused by coupling between vertical movements 
of the pivoted gravity meter beam and the component of horizontal acceleration 
in the direction of the beam (see Figure 2). The centre of mass G of the gravity 
meter beam is at a distance d from the pivot point O. The counter-clockwise 
couple produced by gravity and the horizontal and vertical accelerations #, 2 is 


Md{(g + 2) cos 6+ & sin 6}. (10) 


We will substitute for # and # the expressions given in equation (3). We are 
interested only in the case in which @ is small, in which case the corresponding 
solution for @ can be expressed in an infinite series of sine and cosine terms. We 


yMg + MZ 

Fic. 2.—Cross-coupling in a beam gravity meter. 
will approximate @ by using only the largest sine term but the analysis could be 
extended to include the other terms also. The assumed expression for 6 will 


then be 


+ 6; sin(wet + (11) 


Substituting this equation into equation (10), neglecting 6° and higher powers 
and taking an average over a time long compared with the periods concerned, we 
obtain 


Couple = Md{[g cos 6o(1 — }6;?) 
+ cos sin sin(wgt + 


4 
} 
| 
} 
a vv, 
q \) 
; 
4 
Q 
Q 
4 
q 
| 
0) 
q 
‘ 
(12) 
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If the beam is horizontal in its normal zeroed position and the meter is read 
with the beam close to this position, then 6p is small and the couple becomes 


Md{[g(1 —4 47) 
lg 6,2 

— 2000; cos(y — «) 

+ #0; sin sin(wst 


Of these terms the first is the couple acting on the beam deflected through an 
angle 4 in the absence of wave accelerations and the others are of the second 
order in 9, 6;, with the exception of the last which is zero unless w; = ws, in 
which case it becomes 


(13) 


101 cosy. (14) 


This is the cross-coupling term and is only important when the movements of 
the beam have the same period as the horizontal accelerations. The periodic 
angular movements of the beam are driven by the vertical accelerations, so this 
condition is equivalent to the periods of the vertical and horizontal accelerations 
being equal. In circular motion, probably a fair approximation to actual condi- 
tions at sea, the vertical and horizontal accelerations have equal periods and a 
m/2 phase difference. 

The magnitude of the cross-coupling effect depends on the construction of the 
meter concerned. One method of making a surface ship gravity meter is to greatly 
overdamp it, so that ship motion will not cause the beam to strike the stops. It 
is possible to achieve the same result by using a beam having a large moment 
of inertia compared to its moment of mass about its axis of rotation; i.e. use a 
beam which is nearly balanced about its axis of rotation. A third method is to 
use a stiff spring. The case of an overdamped gravity meter is considered below 
and it is shown that large errors can result when the ship’s motion is circularly 
polarized. Similar consideration of the other two cases would show large cross- 
coupling errors for plane polarized ship motion, which motion appears at points 
displaced both horizontally and vertically from either the pitch or roll axes of the 
ship. 

‘Cecdiine now the case of an overdamped meter and suppose the meter beam 
to be light and the mass to be concentrated, so that the moment of inertia of the 
beam about its pivot point is Md. Let the natural period of the beam be Tp = 27/f 
and the static sensitivity of the meter be S radians/gal. Then the equation of 
motion of the beam acted on by vertical accelerations Z is 


+ BO+CO = dz sin(wst + «) (15) 
ft = Cid. (16) 


B and C are instrumental constants, B is a damping term and C represents the 
elastic restoring couple produced by unit angular displacement of the beam. S 
and C are related by the equation 


S = dC. (17) 
The particular integral of (15) and (16) (the free oscillations are very rapidly 


; 
big. 
| 
i 
a 
| 
{ 
4 
} 
¥ 
3 
a 


Some theoretical considerations in the measurement of gravity at sea 
damped out) is 


6 = + sin(wgt + «—¢) 
d*(f? — 


Comparing (18) with the expression (11) for 6; 


tang = 


= + (20) 
and 
b= (21) 


If the horizontal accelerations lie in the same vertical plane as the meter beam, 
the cross-coupling term becomes 


+ 


The meter must be such that 1 mgal change in gravity can be measured. Suppose 
this needs a movement of 10~4cm at the end of the beam. Then 


10°3 .dS = 10-4cm (23) 
and from (16) and (17) 


dS = 1/f?. (24) 


f? = 10. 


The meter must be damped so that the beam will not hit the case in the presence 
of large vertical accelerations. Put 


amplitude of beam movements at frequency ws Pp 
amplitude of beam movements at zero frequency 
and suppose P = 0-o1 at resonance (ws = f = ¥/10). From (20) and (24) 


and putting f 
W3 = 
= 
Bf/d? = 1000. (26) 


Substituting in expression (22) and noting that ¢ = 7/2 at resonance, the cross- 
coupling term becomes 


cos[e—(7/2)] gal = — cosfe — (=/2)] mgal. (27) 


For the case of circular motion of a meter (€ = 7/2) with 10cm beam length, and 
with #; = 2; = sogal, the error introduced by cross coupling is 


50x 50 


coso = 125 mgal. 


20 


: 
95 
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The phase angle ¢ is 7/2 at resonance, the worst possible value for circular 
motion, but the magnitude of the effect in practice increases with period of the 
vertical accelerations because, owing to the heavy damping, the increase in P 
more than offsets the effect of change in phase angle. Its magnitude (EZ) is shown 
below for several values of m, where f = nwg and for a beam length of 1ocm. 


I 2 


E(mgal) 125 250 500-1243 
¢ (deg) 80° 88° 84° 


The effect is clearly a most important one. It can be reduced by a slow rotation 
of the meter about a vertical axis, which has the effect of varying ¢. It is, of course, 
possible to make a simple computer to calculate the cross-coupling effect and to 
subtract it from the gravity meter reading. Another approach is to keep the 
meter beam fixed in its zero position by a fast-acting feed-back loop which would 
apply magnetic or electrostatic forces to the beam. The average force applied 
to the beam would have to be computed and added to the meter spring setting. 
The purpose of the analysis given above is to show that the effect is present and 
can attain a magnitude of several hundreds of milligals. It therefore cannot be 
ignored. 

The preceding investigations have assumed the wave accelerations and 
stabilization errors to be simple harmonic. It can be extended with very similar 
conclusions to the case where these quantities are considered to be the sum of a 
number of such terms. In the case of stabilization errors the product of tilts and 
horizontal acceleration terms with the same period produce a number of terms 
similar to (7), while products of terms with unequal periods disappear. In a similar 
way cross-coupling terms only arise from the products of components of equal 
period in the horizontal and vertical accelerations. 


3- Theory of gravity meters free to swing in gimbals 
Figure 3 represents a gravity meter free to swing about a gimbal axis O, its 
centre of gravity is at G and the sensing element (the mass) at P. Let OP = /, 
ZOAP = 90°, and the free period of the gimbal system be T» = 27/wo, where 
wo? = g/Lo. The point O is subjected to periodic accelerations %, 2. 
Then the equation of motion of the gimbal system is 


6 + FO + + 2/g) sin + cos = 0. (1) 
The acceleration of P normal to OP (a7) is given by 
ar = 18+ (2) 
and the acceleration of P along PO(a,) is given by 
a, = (3) 
The component of acceleration parallel to GO(a,) is given by 


@g = 
= 16 cos a—#sin 6+ (g+2#)cos 6+ /6sin« 


| 
| 
| 
| 
a 
| 
| 
| 
4 
4 
4, 
i 
4 
: (4) 
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and the component perpendicular to GO is 
Ayn = ay Sin % — ayCos 
= 1@sin«—#cos 0—(g+#)sin cos x. (5) 
¥ may be eliminated from (4) and (5) by substituting from (1) 
dg = 16 cosx+l0sin «+ (g/t? )tan + FO + + sin 0}+(g+%)cos@ (6) 
an= 162 sin x—10 cos (g/wo?)(0 + (7) 
Z 


Fic. 3.—Gravity meter free to swing in gimbals. 


In order to make an instrument to measure gravity at sea, it is desirable to have 
an analogue computer perform the calculations indicated in equations (6) and (7). 
It is therefore very desirable to simplify them. This can be done by making OA 
equal to Lo so that 


1 = g/wo* cos «. 8 
(8) 


ag = (g/two2)(62 + «+6 tan + FO tan 6) +(g+2)(cos 6 + sin tan @) (9) 
an = (g/wo2){6? tan «+ FO}. (10) 
Expanding (9) in powers of 6 
ag = (g/two?)(62 + 00 + 4003 + tan + FOO + 4 F096) + + 62/2 + e464) + +462) (11) 
neglecting terms in 65 and 264. 
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We are interested in obtaining gravity readings over finite time intervals t; 
to t2 (about 10 minutes). Over such time intervals @ can be expressed as a Fourier 
series 


6 = A+ > Businwgt+ > Beycos wy. 


Substituting this expression for @ into equation (11), taking # = o and taking an 
average over a time long compared with the periods concerned, we obtain 


iy = B+ + 00° + (12) 
The correction to be subtracted from the measured acceleration, ay, is given by 
hgh + + + 4262 (13) 


and the second-order term is 
(14) 
The magnitude of the higher-order terms will be evaluated by an approximate 
method in some particular cases in order to estimate the errors caused by neglect- 
ing them. Since we are interested only in the case in which 6 is small, equation (1) 
may be approximated by 


6+ FO + 2920 + (wo2/g)20 + wo2x/g = 0. (15) 
# and 2 will be assumed to have the form 


&/g = eysinwt 


= €,sin wt + €-cos wt. (16) 


As previously noted the solution of (15) for @ small can be expressed as an infinite 
series of sine and cosine terms. This series will be approximated by 


( = A+ B,sin wt + B, cos wt + 5, sin 2wt + 8, cos 2wt (17) 


where B;, B, are first-order terms in €», €s, €- and A, 5, and 8, are second-order 
terms. Fifth- and higher-order terms will be neglected. 
Differentiation of (17) gives 
6 = A+ B,sin wt + B,cos wt + 8; sin 2wt + 6, cos 2wt 
6 = wB, cos wt — wB, sin wt + 2w8, cos 2wt — 2w8, sin 2wt 
6 = —w2B, sin wt — cos wt — sin 2wt — 45, cos 2ut. (18) 


These values are now substituted into (15) and the coefficients of the constant, 
sin wt, cos wt, sin 2wt and cos 2wt terms each equated to zero to give the follow- 
ing equations: 


B 


(wo? — w*)B, — wWFB, + = 0 
(wo? — w?) Be + wF Bs + =o 

(s00? — 400%)3, — oBe 


2 
woerBs 


(wo? — 


= 0. (19) 
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Solving for B; and B, 
wo? [(wo” — w* es + ecw F] 
(tog? — + w2F2 
— we — F} 
+ 


and solving for A, 5s, 5¢ in terms of wo, w, F, Bs, Be, €v, €s and € 


B,; = 


Be = 


— Bey 
2 

[(wo? — 4w*)B,— 2wF 

[(w0? — 420°)? + 4m? F?] 
Wo"€y [(wWo? — 4w*)B, + 2wF Be] 

[(wo?— + F2] 
Substituting (16) and (17) into (13) the following expressions are obtained for the 
higher- than second-order terms, labelled A, B, C. 

= BeBe + (25) 


+ = + (26) 
5;Be 


5; = 


2 2 2 2 


= = - |: (27) 


It should be pointed out that the approximations made in deriving the preceding 
expressions for A, B and C are probably not valid when resonances are approached. 
Methods of avoiding such resonances will be discussed later. 

A, B and C have been evaluated for circular motion in a vertical plane with 
frictionless gimbals. Putting 


(28) 
it is seen that all terms contain a*. They are tabulated below in mgal for various n 
and a = 0-05 their values for other accelerations are equal to the tabulated value 
multiplied by the fourth power of the ratio of the acceleration amplitudes. The 
percentage error increases as the square of this ratio. 


Second-order Percentage error in 
A+B+C term second-order term 


to 


6) 

+6°5 1410 0°46 
@ 

972 

850 

791 

677 

—0°%3 625 


99 
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All the terms are infinite for n = 1 and C is infinite for m = 2. These correspond 
to resonances of the gimbal system when driven by accelerations equal to the 
natural period of the gimbal system or twice this natural period. In the absence 
of damping, @ becomes infinite. However the resonance is sharp and at 1/5 times 
the natural period, the error term A+B+C has already fallen to 5-3 mgal, or 
0°55 per cent of the correction applied. The sum of the higher-order terms does 
not exceed 3-3 mgal if there is no energy of period shorter than 4/5 times the 
natural period, and if the magnitude of the second-order term }g@? is limited to 
600 mgal. The natural period of the gimbal system is likely to be about 1 second. 

The very large values of 0, and of all correction terms, at the resonances of the 
gimbal system can be avoided by introducing damping into the system. A, B and C 
have been evaluated for critically damped gimbals with 


€ = €y = a = 0°05 
€,=0 
F = 2m 
Wo = nw. 


Cc A+B+C 
mgal 
n? = 1 +0'2 —0'4 
-o'l 
4 +0°2 —o'6 
5 +03 —o'6 —0'4 
7 +03 -o'l —o'6 
9 +04 —o'l 
25 +0°5 ° —o'2 
+0'5 ° —o'8 


However the introduction of damping causes difficulties with cross-coupling 
effects, which are absent when the gimbals are frictionless. The acceleration 
perpendicular to the measuring axis of the gravity meter is given by 


An = (g/wo?){6? tan «+ FO}. (10) 


In the case of frictionless gimbals, F = 0, and a, becomes a second-order quantity. 
Substituting from (18) for 6 and neglecting third- and higher-order terms, 


Qn = g(w?/wo")(Bs cos wt — B,sin wt)? tan «. (30) 
For circular motion described by (28) 


= tan « sin2zt 


= (ga?/2(n?—1)) . tana(1 —cos 2wf). (31) 
Movements of the gravity meter beam have the same frequency as the vertical 
accelerations whereas ay given by (31) has twice this frequency. Hence no cross- 


coupling of second-order terms can occur and cross-coupling is a negligible 
factor. 


In the case of a damped gimbal system, 
Qn = (g/wo”){w?(B, cos wt — Be sin wt)? + 


+wF(B, cos wt — wB, sin wt + 2w5, cos 2wt — 2w8, sin 2wt))}. (32) 


(29) 
: 
is 
j x 
t, 
up 
By 
Ag 
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The beam movements are 7/2 out of phase with # for a heavily damped meter so 
that the only term in (32) which is effective in producing a cross-coupling term in 
circular motion is 

(g/t?) By cos wt (33) 

wFB,coswt = — " 
(wo? — 2+ 
For critical damping F = 2we and the component of a, significant in cross- 

coupling is 


4an? 
(34) 


The acceleration effective in cross-coupling for a meter on a stabilized platform 
is a, so the acceleration given by (34) is smaller than, but comparable with the 
case of a stabilized meter when the damping is critical. Damping, therefore, does 
not provide a satisfactory answer to the problem of resonance of the gimbal system. 
The motion of a fair-sized ship does not contain much energy close to the gimbal 
period, but there may be sufficient to cause trouble with an undamped gimbal 
system. LaCoste & Romberg have solved this problem by using frictionless gim- 
bals and allowing the point of support to move relative to the ship in a manner 
which prevents large amplitudes of gimbal swing. The theory of operation of a 
{ravity meter suspended in gimbals as developed above still holds good, only the 
relevant motion is now the motion of the gimbal pivot and not that of the ship. 
As an example, consider a gravity meter suspended in frictionless gimbals, the 
gimbal pivots themselves being supported by a damped pendulum (see Figure 4). 
xX, 


Axes a 


in space port of damped pendulum 


Ss 
“Phttached to ship 


Gimbal axis of 
gravity meter 


Gravity meter 


Fic. 4.—Pivot of gravity meter gimbals suspended from damped pendulum. 


The two pendulums are treated as simple pendulums of masses m, m2 and 
lengths /;, /2. The linear deflections of the bobs relative to axes fixed in space are 
x, and x2 and the support of the damped pendulum is at distance x from a z axis 
fixed in space. The equations of motion of the two pendulums for small oscilla- 
tions are 

+ + ¥2) = (35) 
+ + Rx + — mogxe/le = 0 (36) 
where R is a frictional damping term. 
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&lh = 
= 
Rim, = 2fw 
m2/m, = a. (37) 
Then (35) and (36) may be rewritten, 


+ = 0 (38) 
+ + + = 0. (39) 
Now suppose the ship to be moving horizontally with simple harmonic motion 
so that, 
x = Acoswt (40) 
and let the motions of the upper damped pendulum and the gravity meter be 
x, = A, coswt+ B,sinwt (41) 
and 
x2 = Agcos wt + Besinwt. (42) 
Expressions for #, ¥1, ¥;, and #2 may now be derived in terms of Aj, Bj, Ae, 


Be, sin wt and cos wt and substituted in (38) and (39). Equating the coefficients, 


of sin wt and cos wt to zero in each equation the following equations for Aj, 
B,, Az, Be are obtained: 


— w*A) + (we? = (43) 
+ — w?)Be = 0 (44) 
— w?) A; + 2fww,B, — aw2*A2 = w*A (45) 
— 2fww, A; + (w;? — w?)B) — aw2*Be = o. (46) 
At the principal resonance of the gimbal system 
w= We (47) 
and equations (43) through (46) may be solved to give 
A;(w2) =-A 
B,(w2) = 
A2(w2) = — 


B2(we2) 


The motion of the gravity meter remains of reasonable amplitude and the 
damped pendulum supporting the gimbals swings so as to keep the gimbal axis 
stationary in space. The swinging of the gravity meter is equivalent to a free 
oscillation of steady amplitude and the correction terms given by (13) are valid for 
this motion. The problem of excessively large swinging of the meter in its gimbals 
at the resonance frequency has been removed. 
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4- Conclusion 


These investigations have shown that gravity observations made at sea on a 
stabilized platform can be greatly in error if there is a coherent phase relationship 
between levelling errors in the platform and the horizontal accelerations. Periodic 
stabilization errors need to be reduced below 8 seconds of arc in order to be sure of 
1 mgal accuracy in the presence of horizontal accelerations of 50 gal amplitude. 

Horizontal accelerations of a spring gravity meter produce a couple on the 
beam once the beam is deflected from its zero position. A persistent relation, 
between the phases of the vertical deflections of the beam produced by the vertical 
wave accelerations and the horizontal accelerations of the meter, can easily cause 
errors of several hundreds of mgal. 

Gravity measured by a meter free to swing in gimbals can be adequately 
corrected by a single second-order term up to accelerations somewhat in excess of 
sogal, provided that the gimbal system does not swing with unduly large ampli- 
tude due to resonance with wave accelerations. For appreciably larger accelera- 
tions it is necessary to include higher-order correction terms. The problem of 
very large oscillations of the gravity meter’s gimbal system at its resonant fre- 
quencies may be solved by proper suspension of the pivot point of the gimbals. 


Lacoste and Romberg Company (L.J.B.L) Institute of Geophysics, (J.C.H.) 
Austin, University of California, 
Texas, Los Angles, Cal., 
U.S.A. U.S.A. 
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Analysis of Gravitational and Geometric Aspects of 
Geodetic Utilization of Satellites 


W. M. Kaula 
(Received 1960 September 22)* 


Summary 
Expressions are derived for the first-order effects of any period of 
.any term, U7", of the gravitational potential on the orbital elements, 
plus the second-order effects arising from the interaction of U," with 
U8, the oblateness. The order of magnitude of some daily and semi- 
daily variations is estimated to be + 100m from statistical data on the 
gravitational harmonics. 

A general geometric and statistical treatment of all types of observa- 
tions is developed, with the purposes of obtaining rigorous evaluations 
of orbital and observational schemes, and optimum solutions for geo- 

detic positions, gravitational harmonic coefficients, and orbital 
elements. 


1. Introduction 


This analysis seeks to place the gravitational and geometric variables involved 
in the geodetic satellite problem in the framework of generalized least-squares 
adjustment of geodetic data (Kaula 1961). In generalized least-squares adjustment 
taking into account correlation (Arley & Buch 1950, Brown 1955), condition and 
observation equations are expressed in matrix form as 


CY+MZ =F (1) 
and the solution is made subject to the condition 
Y7W~1Y = minimum. (2) 


C and M are block matrices comprising the coefficients of the equations. F is a 
column matrix of the residuals of the equations. Z is a column matrix comprising 
corrections to parameters, and Y comprises corrections to observations. The essen- 
tial difference between parameters and observations is that only the latter are 
quantities for which there can be given as starting values in the adjustment esti- 
mates with meaningful variances and covariances, which are the content of the 
covariance matrix W. 

The gravitational and geometric variables entering the geodetic satellite prob- 
lem are (a) the constants of integration of the orbit; (b) the gravitation field 
quantities: the Earth’s mass, the oblateness, and a number of harmonic coeffi- 
cients or area mean anomalies to express the anomalous variations; (c) the posi- 
tions of geodetic stations; and (d) observations of the satellite from the stations. 

* Received in original form 1960 August 16. 
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The relationships between (a) the orbital constants and (b) the gravitational 
coefficients will be first derived. Next, these variables will be combined with (c) 
the station (or datum) positions and (d) the observations in equations expressible 
in the form of equation (1). The final phase is estimation of variances and co- 
variances of all these quantities to construct the covariance matrix W in (2). 


2. Gravitational perturbations of the orbit 

Several theories of satellite orbits have been published taking into account 
the oblateness to at least secular and long-period second-order effects (Garfinkel 
1959; Kozai 1959; Brouwer 1959; Musen 1959). Considering the anticipated 
magnitude, in comparison to the oblateness, of the additional small independent 
terms U” of the gravitational field, it is now therefore of interest to derive their 
effects on a close satellite orbit. Groves (1960) has published a derivation of the 
first-order effects of any harmonic term on an orbit, and Musen (1960) has incor- 
porated the expression for any harmonic in the Hansen theory of satellite motion. 
The analytical development herein is felt to be appreciably simpler, and is directed 
primarily to obtaining expressions to be in turn used in deriving the coefficients 
of differential corrections: i.e. elements of the C or M matrices in (1). ‘The theory 
follows most closely that of Kozai (1959). 


2.1 Trigonometric identities 
The following expressions will be used repeatedly: 


cosmx = = A(cosx+jsin x)” 


m 
8-0 


where #& denotes the real part, and j = (—1)! 


sinmx = = Bi —j(cosx+jsinx)”} 


m 
R > sin’x 
8-0 


J arr b 
sin*x cos?x = [- — + 
2 2 


24 20 
d-0 


x {cos(a + b—2c—2d)x +j sin(a + b—2c—2d)x} 


“4g 
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q 
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cosx, 6 odd 


For even, non-zero terms are: cos(a+)x, cos(a + b—2)x, ..., 
1, beven 


odd 


For a odd, non-zero terms are: sin(a + 6)x, sin(a + 6 — 2)x, ..., 


sinx, 6 even 
cos A cos B = cos( A + B) +} cos(.4 — B) \ 
sin Asin B = —}cos( A + B) +} cos( A — B) 
sin Acos B )sin( A+ B) +) sin( 4 B) 
~ }sin( 4 B). 


(6) 


cos Asin A+ B) 


2.2 Transformation of harmom term to orbet-referred (oor dimates 
The general expression for the sphereal harrmmonm term of degree / order 
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CY+MZ F (1) 
and the solution is made subject to the condition 
minimum. (2) 
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tions of geodetic stations; and (d) observations of the satellite from the stations. 
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From the definition of Legendre Associated Functions, 


= —(sin?d —1)! 


even 
m) odd 
comm 


applying equation (10) 
Apply (3) and (4) to (8), and substitute (9) therein 


coe mA « \ 


me 


2.1 Trigonometric identities 
The following expressions will be used repeatedlyy 


cosmx = Remit = A(cosx+jsinx)” 


m 
& ji cos™~&x sin*x (3) 


where # denotes the real part, and j = (—1)! 
sinmx = B{—jemi2} = B{—j(cosx+jsinx)"} 


> ) cos” sin*x (4) 


sin®xcos?s = [- 
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cosx, 6 odd 
For a even, non-zero terms are: cos(a+)x, cos(a+b—2)x, ..., 
1, beven 


sin 2x, b odd 


For a odd, non-zero terms are: sin(a + 5)x, sin(a + b—2)sx, ..., | 
sinx, 5 even 


cos A cos B = $cos(A + B)+4.cos(A— B) 
sin Asin B = 

sin A cos B = }sin(A+ B)+}sin(A— B) 
cos Asin B = }sin(A+B)—}sin(A—B). 


(6) 


2.2 Transformation of harmonic term to orbit-referred coordinates 
The general expression for the spherical harmonic term of degree /, order m is: 


Up = sin $){Ap" cos mA + Bp sin ma}. (7) 


Substitute «—@ = A where « is right ascension and @ is GST. 
cos’mA ={cos m(«—Q) cos m(Q — 6) — sin m(a— Q) sin m(Q —@) 
8) 
sin mA = sin m(a—Q) cos m(Q—@)+ cos m(a—Q) sin m(Q — 6) 


Fic. 1.—Orbit-equator—meridian triangle. 


In the triangle, orbit : equator : satellite meridian, 
cos(a— 2) = 
(9) 


sin ucosi 


sin(a—Q) = 


sing = sinisinu. 


Sie 
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Weveu = 
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From the definition of Legendre Associated Functions, 


cos md, di+m 


(—m)/2, (@m) even 


\(l—m—1)/2, (l—m) odd nl—m—2ts ei nl—m— 


applying equation (10). 
Apply (3) and (4) to (8), and substitute (9) therein: 


cos = )i cos m( + (12) 


sinm’ = > jie m(Q— 6) —j cos m(Q—6@)}. (13) 


Substitute (11), (12), (13) in (7), cancelling out the cos™¢’s, and combining powers 
of sinu 


(l—m)j2, even 
(—m—1)/2, (—m) odd 
I 2l—2t)\(—1 
(al—2t)\(— 1} 


x 


x —jB™) cos m(Q— 6) + (BM +jAM) sin m(Q_—6)} x 


x > ) cos™~*u 


Apply (5) to (14), with a = 1—m—2t+s,b = m-s: 
I (2/—2t)(—1) 


sint-m-2tj x 
2t)!t\(1—t)! 


xB (AP cos m(Q—6) + (Bf sin m(Q—6)} x 


SS 


x (— 1)*{cos(/— 2t — — 2d)u sin(/— 2t — 2c — 2d)u}. (15) 


Since from (5), (12), (13), all parts of (15) are real, the products are all real, and 
any term which has an odd power of j as a coefficient can be dropped because 
there will be another term of opposite sign cancelling it out. 
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Applying equation (6) to the terms in Q, @, u in (15), and taking further into 
account the factor j/-™-2¢+28 between the s and c summations, we get finally as 
the non-zero terms in these arguments: 


for 1—m even: + 

+ Bf (16) 
and for /—m odd: +A)" —6)}— 

— By” cos{(/— 2t — 2¢— 2d)u+ m(Q—6)}. (17) 


When /— mm is odd, then, the factor should appear in the summation, 
or which is in turn multiplied by (— 1)!-™~2¢+* and (—1)*: 
oad 
—m—1)/2 
1—2t)! 
U? = —— sint-™-2tj x 


at)!t\(/—t)! 


m —(l—m)/2, lm even 


gi-2t 


l—-m-2t+s m-—s 


Br even 


A? 


sin{(/— (18) 


odd 


It is desirable to transform (18) so that terms of the same argument 
cos 
cin 8) 


are collected together, and, insofar as practicable, factors depending on the same 
indices are combined. 

Substituting p = t+c+d in the argument necessitates in turn eliminating one 
index from the factors; putting d = p—zt—c seems most convenient. The elimina- 
tion of the d summation must be balanced by a summation with respect to p from 
o to I, which is placed before the t summation so that the arguments will be 
combined. The limits of the d summation places limits on the possible values of 
c:0 < p—t—c < m-—s, so that the lower limit of c becomes c > p—t—m-+s for 
p—t > m-s, and the upper limit of c becomes ¢ < p—t for p—t < 1—m—2t+s. 
Note also that ¢ < p. 

To combine coefficients incorporating the same indices as much as possible, 
shift the (/—m-—2t)! factor in the denominator of (18) to the right of the 


aa 
Ue 

i ae = 
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summation, combine the 2! and 2!-2¢ just after the ¢ summation, and combine the 
two (—1) factors after the c summation. The expression for U” then becomes: 


Br even 
+ po sin((I- (19) 
where 
(2l—2t)! 
l—m—2t+s m—s 
l-m 
m+t > prs, 


o<t<¢i—,p > —, l—meven; 
2 2 


m+t < p—t—m+s 
m 
| spe , odd 
2 
l—-m—2t+s,m+p+t>I1+s 
m+p+t<l+s 


2.3 Series development and integration 

The problem, as in any general perturbation theory, is to find the most efficient 
method of transforming the perturbing function from its expression in (14) or 
(19), differentiating it with respect to orbital parameters, and placing it in the 
equations of motion for integration with respect to the independent variable. 

For machine computation of orbits, the most effective procedure is probably 
to evaluate the coefficients in equation (14) numerically, substitute for 1/r, 
cos(w+v) and sin(w+v) their series expressions in terms of a, e and the inde- 
pendent variable, and generate the series development of U" by the formulae for 
the multiplication of numerical Fourier series (Herget 1948, p. 123). The type of 
orbital elements used and the independent variable (time, or intermediate orbit 
mean, true, or eccentric anomaly, or longitude) would most logically be those 
used for the second- or higher-order theory taking into account the oblateness. 

Conceptually simplest is to use time as the independent variable and to define 
the intermediate orbit by Keplerian elements. It is necessary to transform 
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For the long-period terms, /—2p+q = 0, integrate equation (21) with 
respect to the mean anomaly by using dM = r? dv/a*(1—e*)!, substituting 
[a(1 —e)/(1 +e cos for binomially expanding (1+ cos and 
applying equations (5) and (6): 


2n 
2p)(w +0) +m(Q—6)} dM 
0 


i-1 


Qn 
0 b=0 
d-0 


+ 2p 0+ 2p + dv. (22) 
cos. 


The condition for long-period variation is (/—2p)+(b—2d) = 0. Since these 
are symmetrically placed coefficients in the binomial expansion, combine the two 
terms and substitute (2d+/—2p) for 6 in (22). There is thus obtained: 


p’-1 


in which p’ = p,__ for p < //2 


(23) 


Short-period terms, /—2p+q # 0, are more inconvenient; we must use 
Hansen’s X?'™ (Tisserand 1889, p. 256), where 


H n= —(l+1), m= (l—2p), i = (l—2p+q). 
= 
where 
B 
2p’ —2l\(—1)' ((/—2p'+q')e\" 
> ( r! 2p (25) 
h=k+q, h=k, q<o. 
(—2p'\ 1 
2 26) 


n=k,q n=k-q,q <0. 
=qforp < lz; p' =I-p,q = —qforp > 
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Note that Gipg(e) = and that Gipg is o(el¢!). 
Combine equations (19) and (21): 


l 
Up = > Fimsli) > M, 8) (27) 
p=-0 
in which: 
A™ even 
Soe =| 
— odd 
+ sin{(I— 2p)eo+(1— 2p +q)M+m(Q—6)}. (28) 
A} i—m odd 


Denote the derivative of Simpg with respect to the argument by Save and denote 
the integrals of Simpg and Siave with respect to time, assuming w, M, Q and @ to 
change secularly, by Simpg and S,,,5, respectively; so that, for example, 


Simpa 


Simpa = (1— 2p +q)n+m(Q_—6) (29) 
In the case / even, p = l/2,qg = 0, m = 0, 
= = % = Ap(t—to). (30) 
’ Then 
= M, Q, 8). (31) 


Differentiate (31) with respect to the orbital elements and put the derivatives 
in the Lagrangian equations of motion (Smart 1953, p. 69) to obtain the time rate 


where the perturbation of the mean anomaly 


é of change of the elements: 

F imp ClpaSimpa 

mpd nal¥84/(1— sini | 

Simpa [V(t cot i 

. FimpSt 


t 
M* = { ndt—nt+o, 

0 
o is the mean anomaly at epoch; 


dG 
ine = from (23) or (24). 


24 al 
imp ~ a rom (20); 
= and (33) 
% 


112 W. M. Kaula 


Continuing: 
Fim S; —e2 | 


To obtain the expressions for the integrated changes in the elements, substitute 
Simpq for Simpg, and Sin nq for Simyg in (32) and (34); e.g. 


FimpGipaSimpa (35) 
35 


AQimpg 


Note from (30) and (34) that there cannot exist, to the first order, any secular 
changes in i, a and e; further that, since /—2p+gq = 0 is the condition for long- 
period changes, there cannot be any long-period variations in a. 


2.4. Second-order effects 


The following changes in the results obtained by integrating (32) and (34) 
should be considered as second-order effects of the potential term Uimpg (31): 
(a) changes in the orbital elements due to Ujmpg itself, which in turn affect the 
change obtained for other elements, and (b) changes in the orbital elements due 
to other terms of the perturbing function, which in turn affect the change in 
elements obtained as due to Uimpg. Mathematically, define the perturbation of 
element E;, AEjmpg (cf. Kozai 1959, p. 374), as: 


aE. 
AEumpa = Ai Eumpet+ > dt 
j 


(36) 
= Ai Egmpqt+ A2Eampe- 


where AE; denotes a perturbation obtained by integration of (32), (34) in accor- 
dance with (29) or (30). A change of type (a) is expressed by the first term under 
the integral, and a change or type (b) by the second term. 

The factors @E;/@E;, AyE; are due to terms other than Uimpg; we consider in 
turn those due to U?—i.e. occurring in unperturbed motion, and those due to 
U8, the oblateness. 

The only element which changes in unperturbed motion is the mean anomaly 
M, for which M = n = (Aja-*)'. Since a, by (34), is subject only to short-period 
variations, this second-order variation in M need only be taken into account when 


: 
: 
| 
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the short-period terms are computed. Applying (29): 


on n S; 


3FimpGipaSimpa 
» l-2p+q#o (37) 


= 0, l-2p+q =0. 


(37) is appreciable compared to the second part of the integral of M* by (32), 
and so should be taken into account. From the second part of the AoE jimpq integral 
all results from U? are zero, since the only non-zero A;Ejo9 = M— Mp is secular, 
and taken into account in the integration (29) by the m in the denominator. 

To consider the interaction with the oblateness U}, write the E; in the 
general form: 


Eampa 


AP 

For E;,,, use the indices 7,s in place of p,q respectively. Differentiate Eumpg 

E;,,-, With respect to the elements and integrate with respect to time to obtain 
AcEiimpg- For long-period terms, (/—2p+q) = 6 = +(2—2r+s); abbreviating 

= c, we obtain for AoEjimpg: 


OE 2Eumpa Ag Ar 


2ors 


(38) 


in which kitmpg = 1 or (1—2p) or (l—2p+q) or m or —(1+3)/a, indicates 
may or may not be differentiated with respect to e or 7, and kj,,,, and yj), 
have similar definitions. This integral representation of AcBunce does not apply 
to three cases: (1) The second term does not apply to the secular effect of 


U3: +0] = 0; 


it is accounted for in the denominator of (29). (2) The first term does not apply to 
the secular effects of other even degree zonal harmonics 


U?: = 0; 


> 
f 
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| 
i 
aa 
7 
4 
= 
‘ 


114 W. M. Kaula 


they also are accounted for by the denominator of (29). (3) The entire integral does 
not apply to the interaction between U} and terms U), / even, for which 
(/—2p) = +(2-—2r) # o and g = +s; these terms can be included in (29), but 
are very small, since Goo-2 = Gaze = 0 by (23). For the non-secular AoEjmmpg 
compare to A)E;impq in similar notation: 


(Ar Pumpa 


AiEumpe 


jsin Q-6 
\ 


Compare the amplitudes of A; and Ag, disregarding the k’s and y’s: 


i 


a-6 


2p)—+b+m—— 
n n n n 


A2Eumpa A} 
AiEumpa 


ap)2+b+ 


A® 00728 x 6-377/2 x 398-61/2 


= 


so the second-order effects should be negligible for m # 0. For b # 0, m = 0, 
terms arising from mp = 200, 201, 202 need be considered. The b # 0, m= 0 


| 
A® I I 
c——m- 
Ag I I 
nas! 
c——-m- (2-—2r)— 
n n n 
AY I I 
n 
n n 
Ao 1—2p) 
nad) w 
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terms will all be multiplied by at least the first power of the eccentricity, since 
otherwise they fall under (3) of the inapplicable cases. For 6 = 0, m = 0, only 
terms arising from /mp = 201 need be considered, since G2o-2 = Goze = 0; 
furthermore, the second term in the integral for 6 = 0, m = o falls under (1) of 
the inapplicable cases. 

Kozai (1959) gives expressions including the second-order effects for /m = 30, 40. 
The expressions for b = o and any /m are given here; since @£2919/@Q, 2 j2010/ dw, 
@Ej2010/2M are all zero, and Aaimpg has only short-period variations, only Q, w 
and M are affected. 


= 


x {—(l—2p)5 cos i+m} 


AS 


Acwimpep-) = 


| 3 


n2al+8{(1— + m(Q — 6)}4(1 — e2)? x 


x {(l— 2p)3(5 cos*i — 1)—6m cos 


The expressions for 6 4 o are generally smaller in effect, and appreciably 
more complicated if developed following equation (38); hence, it would probably 
be preferable to develop them by transformations of canonical variables. The 
magnitude of the b # o effects is such that they would be considered as “third 
order” in the theories of Kozai (1959) and Brouwer (1959). 

Note that all the second-order effects for m = o would become larger and 
require special treatment near the critical inclination i = 63°-4, where @ = o. 


2.5 Numerical estimates 

The term of largest effect in a daily period is probably Uj, since U! and U} 
cannot exist, while the condition /—2p+q = © indicates that the daily periods 
arising from U} and Uj} are multiplied by the eccentricity. If / = 4, ¢ = 0, 
then /—2p+q = 0 gives p = 2. From (31) and (27): 


I 


x {— Bh cos(Q—6)+ Aj sin(Q—6)}. 


x 


Fy2 = 
From (23): 14+3e2 


= 


q 

§ 
{ 
: 
: 

ye 
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Apply (32) and (34): 


15(1 + §e*) 
= 
x {— Bl cos(Q— 0) + Alsin(Q—6)} 


(7sin®icosi “2 sinicosi+ coti){~ 8) + 


(7sinei— 
4 sini 


4120 = 


Mio = x 


— 15¢2/2) 
x {— Bh cos(Q — 6) + Aj sin(Q — 4)} 
1+ $e? 4 
x {B}sin(Q—6)+ A} cos(Q—6)} 
€4120 = ©, 44120 = 0. 
For the Vanguard I satellite, a = 8-67 x 108m, e = 0-19, i = 0°598 = 34°-2, 
n = 0-782, & = +0°000889, 2 = —0-000608 (all rad/ksec) (O’Keefe & others 
1959, p. 248). Using these values with 6 = 0-0728, integrate (41) to obtain: 
AQgi2990 = — 22-5 x 10-®{B} sin(Q — 6) + Aj cos(Q —6)} 
= +61-7 x 10-6{B} sin(Q — + Al cos(Q —6)} 
= +1:1 x sin(Q — + Aj cos(Q—4)} 
Aisioo = +7-9x sin(Q—6)— cos(Q—6)} 
Proposed specifications for a geodetic satellite are a = 7-75 x 108m, e = 0°031, 
i = 78°, m = 0-925 for which Q219 = —0-000211. The resulting perturbations: 
AQ4120 = + 16-3 x sin(Q ~ 4) + A} cos(Q — 6)} 
Aw4i20 = —8-0 x sin(Q — + A}. cos(Q—6)} 
AMiio9 = 0-0 x sin(Q — 4) + Aj cos(Q (43) 
= — 4:6 x 10-6{A1 sin(Q—6)— B} cos(Q— 6)} 
To estimate the magnitudes of the A”, B” utilize the degree variances 


0; {Ag} given in Kaula (1959), Pp. 2413 (use of the estimated coefficients Anm, Bum 
in the same paper will give effects which, on the average, will be too small, for 
statistical reasons). Multiplying by +/{K(2/+ 1)(/—m)!/(1+ m)!}, where K = 1 for 
m = 0, K = 2 for m # 0, to convert from normalized to conventional harmonics, 
and by R'+2x 10-5/(1—1) to convert from milligals to (Mm)!*9/(ksec)?, for 
estimated magnitude of coefficients apply (44): 


\ 
(41) 


[5 
dtlas 4 
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Table 1 


Magnitudes of harmonic coefficients, (Mm)!*8/ksec)? 


+0020 +0859 
0°243 0-382 
0°045 0°024 
+0'018 +0008 


Table 2 
Magnitudes of gravitational effects on Vanguard I orbit 


42 
m 


F115 +30 
+4 

F105 +289 
+92 
+74 


+H 


Table 3 
Magnitudes of gravitational effects on proposed geodetic orbit 


AQ Aw i Aa 
m m m 


+27. 4243 F250 
+4 +4 +5 

o F814 
+69 +34 ° 
+34 +63 ° 
+9 +42 


+ 
N 


Tables 2 and 3 set forth the results of applying the values in Table 1 
(plus of Aj} = + 5-83) to the two orbits, integrating equations (32) and (34). 
The perturbations in metres were obtained by multiplying the perturbations in 
radians by the semi-major axis. Note that: (a) the effect decreases, usually, with 
increase in degree, /; (b) the effect usually decreases with increase in order, m, 
but for some high inclinations it increases; (c) for odd degree harmonics the effect 
increases with increase in eccentricity; (d) the effects Aw and AM* largely cancel 
each other out. 

The combination / = 3, p = 1 or 2, g = 0, indicates that there can be appre- 
ciable perturbations of frequency n. To obtain a statistical estimate of perturba- 
tions of frequency 2 or higher, the rotation of the Earth can be ignored and the 
perturbations assumed to be sinusoidal accelerations of argument 2M, 3M, etc. 
For the amplitudes of these sinusoidal accelerations, use again the degree variances 


o7{Ag} and convert the auto-covariance function from Legendre polynomials to 
cosines: 


cove(dg) = {Ag} conf. (45) 
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The radial, transverse and tangential accelerations thus obtained were applied 
to a circular orbit of semi-major axis 8-ox10%m, using formulas given by 
Routh (1905, p. 256) to obtain the values in Table 4. 


Table 4 
Very short period gravitational effects on orbit 


f i transverse tangential 
m m 

2 +9°0 +29°0 

3 9°0 

4 1°5 

5 

6 


It is thus safe to ignore perturbations of frequency higher than 4. 


3- Adjustment of observations 
General treatments of the geometry for all types of satellite observations have 
been set forth by Veis (1960) and by Henriksen (1960). Some of the same results 


will be derived in this section, as part of the unified solution of the datum, orbit, 
and gravity harmonic coefficient problems. 


3-1 Coordinate systems and conventions 


Denote origins by suffixes: geocentric: G, or none at all; topocentric: 7; 
satellite-centred: S. 


A positive rotation is a rotation which appears counterclockwise to an observer 


looking from the positive end of an axis towards the origin. Position and velocity 
vectors are written with the x coordinate on the first line, y coordinate on the 
second line, and z coordinate on the third line. Rotation matrices are to be sym- 
bolized with the axis of rotation as a sub-index and the angle as an argument. For 
example, positive (counterclockwise) and negative (clockwise) rotations of angle @ 
about the Z-axis are, respectively: 


+cos6 +sin? o +cos6 —sin@ o 
R;(6) = }—sin@ +cos@ o|, Rs(—@)=)+sin@ +cos@ o| = R7(8). 
° ° I ° ° I (46) 


For the differential of a rotation matrix, use the notation: 


dR;(6) = R;(¢) d® 
—sin6d@ +cos@d0 o —sin@ +cos@ 0\/d0 o 
—cos0d@ -sin@d@ o\| -—sin@ o\/o dé 
° 


° ° ° ° o dé 


Since differentiation multiplies every element by dé, the diagonal matrix d@ may 
be shifted to any place in a product of several matrices: in particular, if the product 
forms a vector, to the end of the row as a single element. 


OAC 
J 
wig 
= 
| 
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Pertinent coordinate systems and the rotation matrices between them are 
given in Table 5. Vector designations are only given for those coordinate systems 
to be used explicitly. An alternate rotation matrix notation is Ri, meaning the rota- 
tion matrix transforming coordinates from system j to system i: X; = RyXy. 
For example, in the table Re,; = Ra(—w); Rre = Rai = Raf — Q)Ri( —1)Ra( — 
= RJ(Q)R7(i)RZ(w). In Table 5 all symbols for angles have usual definitions: 
w: perigee argument; 7: inclination; 2: longitude of node; 6: GST; A: longitude 
from G; ¢: latitude; A: azimuth, clockwise from north; z: zenith distance; q: 
parallactic angle; 5: declination; «: right ascension. The sequence of rotations 
from system No. 6 to No. 13 can be considered as an excursion around the astro- 


nomic triangle by the Z-axis, sweeping through each side and turning around each 
vertex in sequence. 


rcosv “ 
Q= sin : Orbit referred; X= v): geodetic; 


z w 


E Xa (48) 
N _ }: Station, or local; A = { yq): Altazimuth instrument; 
H 


p 


XB 
| Station-object vector; B = be: Camera axis. 
ZB 


p 


There is no necessary connection between the geodetic coordinates u, v, w 
and the station vertical direction ¢, A; the latter may be, for example, the ‘‘electri- 
cal” vertical of a radio interferometric system, at an arbitrary direction with 
respect to either the gravity vertical or the ellipsoidal normal. In the Av and 
Py vectors, p is the distance from the station to the object. In the B vector, xg 
and yg are the coordinates measured with respect to the plate centre, with the 
Y-axis determined by the line from the centre to the image of the north pole, multi- 
plied by z2/f, where f is the focal length. If the object direction ap, dp is given 
but not the camera axis direction «g, 5, then the data must be treated as if the 
B and P vectors were identical—of effect only on the weighting, and hence negligible 
in importance if the angle between B and P is a few degrees or less. 


3.2 Perturbation of satellite position 

The osculating elements of the orbit are functions of the intermediate orbit 
elements (or any other set of 6 parameters) and the gravity harmonic coefficients 
(plus drag parameters, luni-solar effects, etc.). Hence a small correction to the 
osculating elements, dE, can be expressed as a linear function ot corrections to the 
intermediate orbit elements, dE, corrections to the gravity coefficients, dA”, 
dB? : 


dE = (49) 
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In (49), J is a unit matrix plus nine elements proportional to At = t—to: 
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em 
“tin bey 
jw 
Bin be. 
@M* @M* 
Gin Bag 
sin 2iAt cosiAt 6A8e cos iAt 
|  1§A8sin 2iAt 21A9(1—5 cos*i)At 3A8e(1 — 5 cos®i)At 
| sin 2i(1 + €2)At _ {3 21A9(3 sin*i—2)(1 3(3 sin?i — 2) A8e(5 + 


2a 


4na*(1 —e2)3/2 8na®(1 e?)3/2 e2)5/2 
The elements of J proportional to A} attain a magnitude of order unity about 
ten days from epoch. The elements of Czg are the @E/@A™, 0E/@B™ obtained 
by integrating equations (32) and (34) and then summing over p and g. For 
example: 
éQ I 
nal*34/(1 sini 


Finp > Gipa* 
p~0 


—gin)-m even 


(l—2p)a + 2p q)n+m(Q—6) 


To obtain the correction to the position vector of the satellite, Q, as defined 
by equation (48), obtain the eccentric anomaly E from M by iteration of Kepler’s 
equation and use (Herget 1948, p. 82): 


(50) 


r cos a(cos E —e) 
dQ = rsinv\= dj ay/(1—e*)sinE 
° 
1—ecosE 1—ecosE 
= ay/(1 —e?) cos E cos E e 
/(1—¢?)sin E, -_*_) 
\ ° ° ° ) 
da 
{ae 
de 


In Q, 2g has no significance, since shifts out of the orbital plane are accounted 


a 
q 
| 
j 
> 
| 
; 
| 
(51) 
= 
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for by dQ and di; Q has been kept as a three-dimensional vector only so that the 
rotation matrices Rg{—w) and R;(—7) can be written following the convention of 
equation (46). 

To obtain the complete perturbation of the satellite position, write from 
lines 1-4 of Table 5: 


= = RxgQ. (52) 


Differentiate (52) with respect to Q, i, w, and Q in turn, and substitute from (51) 
for dQ to obtain: 


and similarly for me @Rre/éw; Cog is the 3x3 coefficient matrix given 
by (51). 


3-3 Observation equations 

Any type of observation can be expressed as a measurement of one or two of 
the coordinates of one of the vectors in Table 5. To transform from the vector 
to the actual measurements, use an observation matrix N, 1x3 or 2x3. For 
corrections to the measurements, denote the transformation matrix by N’; if N is 
not a function of the vector components, N’ = N. For example, for measure- 
ments a and vector A: 


in which 


a= NgsA ( ) 
da = NiadA. 


For photographic observations, the measurements are the plate coordinates 
xp, Yo and the vector is that of the camera axis, B: 


= b = NosB, (56) 


° 


f 


Nos = 
° 


where f is the camera focal length. 
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Normally the image will be close enough to the camera axis that the third 
column of N,,, can be set equal to zero. The direction of the camera axis a, 5p 
and angle of swing gg about this axis necessary to put the X-axis of the plate 
in the camera meridian are established by the camera calibration. The rotation 
Rs3(q) could be applied to the B vector, rather than Ra(—) to the plate coordinates 
to obtain the b vector, as is implied by equation (56). 

If only the observed right ascension, declination, and time are provided, then 
the camera axis must be assumed to be identical with the station—object vector, 
and the observed {xp, yp}ovs. = {0, 0}. 

If the camera is on the ground and the object is a satellite in orbit, the observa- 
tion equation referred to the camera axis will be: 


bons. + dbons. = Beomp. + dDcomp. 
= 4B 
= 
= —9)U0] +N, ,Rex[dXo —Ra( — dUp]. 


where Rgx = from Table 5, @ is Greenwich 
sidereal time, and Up is the assumed or estimated geodetic position of the camera. 

In forming dbcomp. (or dPcomp.), Rex is not differentiated with respect to the 
observed «ag, 5g (or «p, 5p) which it contains, since these directions are here being 
used to define directions of coordinate axes which are essentially arbitrary and 
remain fixed throughout the computation. For the observation equation written 
in terms of the P vector the «p, 5p computed from X7 could be used, in which case 


(60) 


{xp, yP}comp. = {0,0}, {xp,yP}ons. {0,0}, 


but there would be no point in doing this extra computation. 

Equation (60) differs from the general photogrammetric problem in that (60) 
assumes that the camera orientation is fixed. This is a valid assumption if the 
orientation was determined by calibration with several stars; furthermore, we are 
not interested in camera orientation, and (60) would still apply even if the orienta- 
tion ag, 5g was incorrect, so long as the observations bons. = {xp, yo} precisely 
reflected the relationship of the observed object direction ap, 5p to the arbitrary 
camera axis direction Sz. 

Equation (60) constitutes the observation equation for the simultaneous tech- 
nique, which has been most completely treated by Brown (1958a, 1958b). 

For the more general case where the orbit is utilized, the X¢ in (60) must 
be computed by a rigorous orbit theory using the assumed gravitational field 
expressed by two (Aj = GM, A}) or more (Aj, Bi”) coefficients, and equation 
(53) must be substituted for dX g so that the observation equations become: 


bobs. + dBovs. = Beomp. | —Rg(—6)dUo|. (61) 


If the camera is located in the satellite and the object on the ground, equation 
(60) becomes 


bons. + dbons. = NosRaxXs+N,,RaxdXs 
= NosRex[Ra( — 9Up — X dUp— dXc]}. (62) 
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In (62) we are again assuming that the camera orientation has been fixed by 
calibration ; in this case, it would have to be fixed by another camera in the satellite 
taking a simultaneous photograph of the stars in a known relative orientation to 
the camera photographing the ground. 

Placing the camera in the satellite introduces two new difficulties: (1) there 
will be many objects of known geodetic location Up, instead of one; (2) photo- 
graphs taken from different locations X q will include the same objects. These 
characteristics make the problem one of the stereo-triangulation, which seems best 
solved as a separate adjustment by either analytic techniques (Brown 1958a; 
Schmid 1956-1957) or by an image-matching and analogue method designed to 
cope with the large number of images. The output of such an adjustment would 
be a set of satellite positions U g dependent upon the datum Up of the objects of 
known geodetic location. The observation equation thus includes all three com- 
ponents of the position vector and could be written: 


Ug+dU g+dUp = 
dA™ 
dB” 
For theodolite observations, the corresponding vector in equation (48) is A. 


The computed direction is most conveniently obtained from the L vector compo- 
nents: 
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The observation correction matrix is then: 


° 


And the observation equation: 
+ daops, = Acomp, +N, dX ¢—dUpo], (66) 
in which the substitution from (53) may be made for dX g, and 


Rav = Ro( ~)Re( 
from Table 5. ‘The ¢, A used in Ray are the astronomic position, defining the 
coordinate axes with respect to which the z, A were measured by the theodolite. 
Radio interferometric observations such as Minitrack measure direction cosines 


with respect to the axis established by the antennas. If this axis is in the north— 
south direction: 


Ni = o} 


\ Zcomp. H 
@eomp. =| E (64) 
comp. 
tan-1— 
N 
aA = . ( 5) 
” 
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in the east—west direction: 


in azimuth A;: 


(67) 


where p is the distance from the station to the satellite. ‘The observation correction 
matrix is then: 
N’ sin ENcos A; cos A; ENsinA,+N®cos A; 
H(Esin A,+ N cos Aj) 


If interferometric observations are used only in the vicinity of 
EsinA; Ncos A; 
l= + 


p p 
then we can set Njz = Niz. The observation equation: 


lovs. + dlons. = (69) 


where, from Table 5, 


The A), ¢ and A used in these equations should be those obtained from the calibra- 
tion of the interferometric system. 


The geometry of radio interferometric observations has been treated in more 
detail by Kahn (1960). 
For radio ranging measurements, the P vector in (48) would appear to be 


applicable. However, to avoid calculating ap, 5p, use of inertial coordinates is 
more practicable (Henriksen 1960): 


p = 
= ([x— uo cos 6+ vosin 6)? + [x — uo sin vo cos 6]? +[x—wo]*)*. 
From which 


I 
dp = R3(— 6) dUp], 


XT YT 2T 
N,x = (71) 
The observation equation: 
Poos. + dpons. (X7X) + [dX — 6) dUp]. (72) 
Radio interferometric and range observations have been written in equations 
(69) and (72) as occurring instantaneously. However, radio observations are 
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continuous, and the fullest utilization thereof would require several equations 
such as (69) and (72) for each observed passage of a satellite. Before using them 
in an adjustment with other observations, it would be best to combine these 
equations—i.e. to integrate them with respect to time. Using the notation 


tdt = hus, 
te 

etc., rewrite (69) and (72) as: 


(73) 
(Ni udUo 
Bisons. = Pyicomp. + Ral —9))ij dUo. (74) 
Writing out some of the integrals, e.g. 


(NizRz vRs(0)X 


4 [(—sin Aysin A—cos Aysin ¢ cos A)(x cos + y sin 8) + 
+ (sin A; cos A— cos A;sin ¢ sin A) — x sin + y cos + z cos A; cos 
{[*— (wo cos 8+ vpsin + (75) 


leads quickly to the conclusion that the integration should be performed numeri- 
cally for the five minutes or less duration of t;—%;. Evaluation of equation (32) 
with A} = —17°55 (mx 10®)4/(sec x 103)? indicates that significant perturbation 
of the orbit is caused by the oblateness in a time on the order of five minutes, 
while Tables 2, 3, 4 indicate that the effects of other gravitational coefficients 


Fic. 2.—Orbits to be distinguished by range observations. 


should be 15 m or less in five minutes. Since the latter effect is beginning to 
be appreciable, it probably is better to break up the integrals into intervals over 
which dX g can be considered constant—of the order of two minutes. Cutting 
down the interval also is desirable from geometric considerations; for example, 
in Figure 2 an integration of range from t; to tg would not make much distinction 
between orbit A and orbit B, while integration from ¢; to tz or from tz to ts would 
make a great distinction. 
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The numerical integration is probably best performed with x, y, z obtained 
directly from the general programme, rather than going back to the equations of 
motion. 

The final type of observation, Doppler, is also continuous, and hence best 
treated by integrating with respect to time: 


ty 
J Fons. dt = 


ti 
(p5—pi)ovs. + 4(pj—piovs. = + dK 
dX — 9)};— {N,~Ra(— 9)};] dUo. 


(76) 


In this case, it suffices to evaluate equation (72) at the starting and ending times 
of the observation. Again, it is desirable to divide the interval of integration for 
geometric reasons; for example, to distinguish between orbit C and orbit D in 


Figure 3. 


Fic. 3.—Orbits to be distinguished by Doppler observations. 


3-4 Covariance matrix 

Since the relationships are not simple between the different variables entering 
in the observation equations, it is necessary that realistic, and complete, estimates 
be made of their statistical characteristics, and statistical procedures applied in the 
computation, if optimum solutions are to be attained, meaningful comparisons 
made between alternative orbits and observational schemes, and discrepancies 
caused by errors or omissions in physical theory distinguished from those arising 
from inaccuracy or incompleteness of observations. 

Of the observations proper, the simplest problem is probably that of the 
camera observations (equations (60) and (62)) in which the errors are relatively 
uncorrelated between different exposures. The error is principally that of the 
plate location of the satellite image, for which Brown (1958b, pp. 3-4) states + 3u 
is a typical value. The principal cause of correlation between errors would prob- 
ably be timing, and, comparing its effect to that of the plate location error at typical 
focal lengths, such timing error correlation should be taken into account if it 
exceeds 0-5 ms. 

For utilization of satellite-borne camera observations (equation (63)), the 
covariance matrix of the U g resulting from the purely geometrical stereo-triangula- 
tion adjustment is needed for their effective application to the dynamical problem. 
A description of this matrix is given by Brown (1958a, pp. 26-31); usually the 
off-diagonal elements, which express covariance, will be of appreciable effect. 
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To use observations which occur several times in a single satellite passage, such 
as theodolite directions, or which are continuous, such as the radio observations, 
the auto-correlation function, or power spectrum, of their errors is needed. The 
error at any instant can be represented as: 


c= dk = (77) 
Then for the integral with regard to time from f; to fz: 


k=0 


and for the mean square error of the integral with regard to time: 


2a 2 
0 k 


Equation (78) indicates that the mean square error increases with increase 
in the interval of integration; however, the square of the integral increases at a 
faster rate, so that the proportion of error to total magnitude decreases as it should 
with more information. If it is desired to make the variances for different length 
intervals comparable, then the variances as computed by (78) should be divided 
by (t2— 

For the error variance and covariance of instantaneous observations, such as 
single frames of a photo-recording theodolite, from (77): 


var{e} = m{e2} = 
k 


cov{e,te—t1} = mf{eye2} = cos k(t2— 
k 


If several theodolite observations of a passage were summed—the discrete 
equivalent to the integration of the radio observations, equations (73)-(75)—the 
error variance of the sum of n observations at spacing At will be: 


The statistical parameters ay in (78)-(79) would have to be determined by 
harmonic analysis of the discrepancies from camera observations in tracking an 
object with the same significant characteristics as the satellite. In the case of the 
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theodolite, the elevation angle, traverse rate, and illumination are of significance, 
so the test could be accomplished with aircraft. In the case of the radio observa- 
tions, ionospheric refraction is also significant, so that a fully certain test could 
only be carried out with an actual satellite. This ionospheric effect also could 
cause appreciable covariance of observations of different satellite passages or 
from different stations, due to inadequacy of the ionospheric model used in 
reducing the observations. 

The other variables in the observation equations, dUo, dE, dA)", dB", are 
considered as corrections to parameters—i.e. part of the Z matrix in equation (1) 
—if no account is taken of any observational information about the Uo, Eo, AP, 
B* other than that which appears on the left-hand sides of equations (60)~(63), 
(66), (69), (72)-(74) or (76). However, since it is desirable to use as much informa- 
tion as possible (or, in terms of the adjustment, to maximize the degrees of freedom), 
some, or all, of these variables will enter the equations with initial estimates which 
we want to influence the final solution. In this case, they must be considered 
mathematically as observations: as part of the Y matrix in equation (1), which 
means that corresponding elements in the covariance matrix W in equation (2) 
are needed. 

Estimates of Up may come from previous adjustments of the same form, either 
of terrestrial geodetic data (Kaula 1961) or of satellite data, or from conventional 
astronomic position observations. If the Up (and A,", B”) were treated as obser- 
vations in the previous adjustment, then the necessary parts of the new input W; 
matrix are simply lifted from the corresponding rows and columns of the old 
output covariance matrix Wo. If some of them were treated as parameters in the 
previous adjustment, then utilize the rule that the covariance matrix of a linear 


transform Y; = Ag, Y;+ B, is obtained by pre- and post-multiplying the previous 
covariance matrix by the transformation matrix: W, = A,;,W,A™ (Brown 1955, 
p. 15). In the present case, 


| 


(82) 


If the Up is obtained from an astronomic position, then, for an area of average 
geophysical variability, o2{E} = o2{N} = (5-7")? = 30000m?, from Kaula (1959, 
p. 2418), and o?{H} is the sum of the variances of the geoid height and ellipsoid 
radius—something less than 1o000m?. Most Up obtained from astronomic posi- 
tions will not be in average geophysical areas, however, but on islands. The 
best sample for estimating the variance of astronomic positions on islands is 
constituted by 69 islands in the West Indies connected through Hiran trilateration 
to the western hemisphere geodetic net; they yield 


o2{N} = o®{E} = (11°5")? = 124 000 m2, 


Estimates of Aj", Bl” obtained from a previous adjustment would have a 
covariance matrix obtained in the same manner. If the previous adjustment was 
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of terrestrial geodetic data, then A} = GM is a function of two variables— 
= mean radius, R, and the mean gravitational acceleration on the geoid, Go; 


= GoR? (neglecting the terms of order 1/300 due to rotation and oblateness), 
so ae the error variance 


var{e(Ap)} = var{e(R)} + 4R®Go cov{«(R), «(Go)} + R4 var{e(Go)}. (83) 


For lm # 00 or 20, a permissible estimate of A” or B;" is 0, since the actual 
magnitudes of these coefficients are small enough that their effects can be con- 
sidered as linear. Hence, even when we do not want to use previous information 
of A”, B”, they can be treated as observations. The variances in this case will 
then be based on the total variance for the harmonic degree /, as computed by the 
square of equation (44), and the covariances will be zero. To minimize the 
possible effects of statistical irregularities in the o?{Ag}, deliberate overestimates 
can be used, such as the values given in Kaula (1959, p. 2419). 

Two possible cases exist where the intermediate elements Eo might be treated 
as observations: (a) estimated Ep are carried forward from a previous epoch; 
(b) estimated Eo are obtained for the same epoch from observations some of which 
may be included in the adjustment. Case (a) would occur if the adjustment was 
being applied to the same orbit for successive periods in turn, in order that the 
matrices in the computation be kept below a certain maximum, or that small terms 
due to drag or other effects could be neglected. Case (b) would occur if the 
advantage of the more accurate elements from observations over a longer period 
is to be utilized in the adjustment of several observations in a short period, such as 
might be done for datum correction. In all these cases, it is necessary that there 
be full consistency between the estimated Ep and the estimated Uo, Aj", B;” 

For case (a), where the elements Ep at epoch fo are obtained from hase at 
epoch 


Ey = E.1+ 2, l= dt +luni-solar terms + drag effects, etc. (84) 


Assuming E» comprises mean elements, and only secular change due to gravita- 
tional terms are considered, from (30): 


By = Fit = Bat (Bs) 
. 


leven 


where the @E1,29/0A? are obtained from (32). 
Hence for a linear transformation of the same type as (81): 


Wo = ,. 


If allowance is to be made for neglected drag terms, etc., expressed by para- 
meters k;, ke,... then their approximate effect on the orbital elements should be 
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estimated and a corresponding addition made to the W matrix; i.e. 


k 
= Bex ke = BexK, 


AW = (87) 


where Wx will usually be a diagonal matrix with the estimated mean square of 
the k’s as elements. 


In case (2), equation (81), etc. could be applied as soon as the various and co- 
variances of the observations, datum parameters, gravitational coefficients and 
orbital elements resulting from the adjustment over the longer period are known. 


3-5 Adjustment procedure 
To adjust r observations and s parameters to ¢ conditions: 


C Y+M Z=F 


R = YTW;'Y = minimum. 


Set (Kaula 1961) 


T = 
S = 


= 

Z = M;*(Fi-C1Y). 
The output covariance matrix of the observations: 

Wo = Wi-— (TW: T7) 

The covariance matrix of the parameters: 

V = 
The quadratic sum: 

R = S7(TW,T7)"'S. 


ae 

M; F, 7 
; sxr sxe 8x1 
as (89) 

Then 
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Note that it is also of computational convenience to treat all variables as 
observations, to avoid the operations expressed by (89). 

Assume, for example, there are j = 1, 2,...m camera observations from 
h = 1,2,...p datums of g = 1, 2,...k satellites, and gravitational effects up to 
degree L are taken into account, so that there are f = 1, 2,... [((L+1)?—5] gravity 
coefficients (omitting the inadmissible 1st and 2nd degree terms); assume further 
that all variables are treated as observations except orbital elements, which are 
treated as parameters. In this case, there would be 2n observation equations 
involving 22 + 3p + [(L. + 1)? — 5] observations and 6k parameters; to adapt equation 
(61) to the form of (88) set: 
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2nx3p 


(95) 


The other types of observations would be similarly treated. 

As indicated by equations (89), (92) and (93), the output covariance matrices 
Wo and V are not functions of the F matrix, which, as shown by (95), incorporates 
all observational data. Hence (92) and (93) can be used to evaluate different 
observational schemes prior to observation: to estimate the effects of changing 
orbital elements, or the number and location of tracking stations, or the type and 
accuracy of instrumentation, or the frequency of observation. 


4- Conclusions 

The foregoing analysis used with an orbit theory which rigorously takes into 
account the oblateness and luni-solar perturbations, appears to suffice for the 
geometrical and gravitational aspects to obtain geodetic results for a few days’ 
duration from non-peculiar orbits at altitudes in excess of 1oookm. A com- 
plete treatment of observations must also consider the problems of refraction, 
aberration and timing. For lower altitude orbits, consideration of non-linear drag 
effects is needed. For longer durations,: radiation pressure, relativity, tidal varia- 
tions and long-period drag variations must be treated. Any of these effects for which 
an adequate mathematical theory exists can be represented by parameters or 
observations in an adjustment of the type (88)-(94). 

Further special treatment is needed for peculiar orbits which approach the 
conditions of circularity (e = 0), resonance [m/(l—2p+g) ~ n/(6—Q)], or critical 
inclination (a = 0, or sin®i ~ 4). 
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Seiches in a Narrow Lake, Uniformly 
Stratified in Three Layers 


N. S. Heaps 
(Received 1960 September 26) 


Summary 

This paper gives a theoretical analysis of the longitudinal seiches 
(or standing waves) in a narrow thermally stratified lake. A mathe- 
matical model of the lake is postulated consisting of three horizontal 
layers of water, density varying uniformly with depth in each. The 
model gives a good representation of the stratification encountered in 
nature, and leads to a relatively simple analytical solution of the 
hydrodynamical equations. 

The theory is applied to find the periods and modes of ordinary 
and internal seiches in the northern basin of Lake Windermere. The 
results show that in addition to seiche oscillations centred at thermo- 
cline level (predicted by a model of the lake consisting of two hori- 
zontal homogeneous layers of water) there are modes of oscillation 
in which the vertical displacement of water is greatest in the lower 
depths of the lake below the thermocline. It seems likely that the 
presence of such modes could explain the vertical oscillations of 
temperature, of large amplitude and complex wave form, observed in 
those depths. Further investigations are required to test this hypo- 
thesis. 


1. Introduction 


Recent experimental and theoretical studies of water movements in Lake 
Windermere during the summer season of thermal stratification (Mortimer 1952a, 
1952b) have shown that the seiche motions or standing waves derived from simple 
two-layered, or three-layered, theoretical models of the lake (in which the lake is 
supposed to consist of two, or three, horizontal homogeneous layers of water) do 
not adequately describe observed movements of water in the lake, particularly in 
the lower depths below the thermocline. These movements appear to be influenced 
by internal seiches which are not represented in the analysis of the two- and 
three-layered models. Mortimer concluded that a complete interpretation of the 
observed movements must await a theory of internal seiches in a basin in which the 
density of the water varies continuously with depth as in the lake. 

A number of papers have been written dealing with internal progressive waves 
in open continuously stratified water (e.g. Fjeldstad 1933, 1936; Groen 1948), but 
the only well-known theory of internal seiches in a continuously stratified lake 
appears to be that given by Wedderburn (1912). This theory, however, is of 
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limited application, for it can only be expected to give an account of seiche oscilla- 
tions at thermocline level, in view of the assumption of a “boundary of zero 
horizontal motion”’ at this level. The apparent need for a more general theory of 
longitudinal seiches in a narrow lake in which density varies continuously with 
depth, which gives an account of oscillations at all levels in the lake, is met in the 
present paper by considering a lake of rectangular longitudinal section (Figure 1) 
in which density varies linearly with depth in three horizontal layers, being con- 
tinuous at the interfaces between the layers (e.g. see Figures 2 and 4). Vertical 
density gradient, constant within each layer, but differing between one layer and 
another is discontinuous at the interfaces. Experimental observations of tempera- 
ture in Lake Windermere (Mortimer 1952a, 1952b) and in other lakes (Wedder- 
burn 1912, Mortimer 1953) during the summer season show that the assumption 
of three layers with a constant vertical density gradient in each is reasonable, 
i.e. in making the assumption a good theoretical representation of density condi- 
tions in a natural lake is obtained. 

The significance of the present work is that it shows how, by an easy mathe- 
matical analysis, an analytical solution to the problem of internal standing waves 
in a stratified lake may be derived without any serious misrepresentation of the 
vertical density distribution in the lake. An analytical solution, giving periods 
and modes of oscillation in terms of parameters of the density distribution, was 
sought, in preference to a solution obtained using numerical integration (e.g. 
Fjeldstad 1936), in order to obtain an understanding of the mathematical structure 
of the motions. An alternative approach to the problem is to consider a theoreti- 
cal model consisting of more than three homogeneous layers of water of differing 
density (Makkaweev 1936; Proudman 1953, 359). As the number of layers is 
increased the more nearly can the model represent the oscillations in a lake with 
a continuous vertical density profile. The complexity of the equations governing 
the oscillations mounts rapidly, however, as each additional layer is introduced. 

The work in this paper forms part of the solution of the larger problem aimed 
at the complete determination of water movements in a lake subject to wind 
action during the summer months of thermal stratification. As described by 
Mortimer (1952b) the response of a lake to a wind impulse is firstly one of forced 
movement under direct wind stress, followed by one of free oscillation in the form 
of surface and internal standing waves, or seiches. The present work sets out to 
determine periods and modes of free oscillation. This is a necessary preliminary 
to the determination of water displacements and velocities which will, in general, 
be the result of the superposition of several modes of oscillation (e.g. see Section 8 
below). 


2. Assumptions 
Assumptions made in the theoretical analysis are summarized as follows: 
(i) The lake, in equilibrium, consists of three horizontal layers of water, 


density varying linearly with depth in each. The equilibrium depths of 
the layers do not vary with time and are, in general, different. 


(ii) The lake is long and narrow. Its longitudinal section is rectangular and 
the depth is small in comparison with the length. 


(iii) There is no motion of the water in a transverse direction: the motion is 
two-dimensional in a longitudinal section of the lake. 
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(iv) Velocity components and vertical displacements, together with their 
derivatives with respect to distance and time, are small enough for terms 
involving their squares and products to be ignored. 


(v) Vertical acceleration is negligible. 
(vi) Geostrophic effects may be ignored. 


(vii) The viscosity of the water, forces of friction, and turbulent mixing may 
be ignored. 


(viii) The water is incompressible. 


Although vertical displacements are often large in practice, the assumption 
that they are small is necessary if the equations of the analysis are to be success- 
fully linearized, so that the motions reduce to simple harmonic type oscillations. 


3- The equilibrium density distribution 
The length of the lake is /. When the water is in equilibrium, the depth of the 
lake is A and rectangular coordinates (x, y, z) are taken with origin in the surface 
of the water at one end of the rectangular longitudinal section, the x-axis being 
_ horizontal (along the length) and the x-axis vertically downwards (Figure 1). 
In equilibrium, the vertical density distribution of the water in the lake 
(Figure 2) is postulated by 


foro<z<h 


for hy < z < hy+he 


my for < < A) / 


where po is a constant density and a), az and ag are constant lengths such that 
hy/a1, he/az and hg/ag are small in comparison with unity (since density differences 
in a stratified lake are small—of the order of one part per thousand). 

Thus, in each of three horizontal layers of depths /, Ae and hg, the density of 
the water varies linearly with depth (p0/a1, po/a2 and po/az are the density gradients 
in the three layers respectively) and is continuous between one layer and the next. 


4- Basic equations 
In the free oscillatory motion to be considered let u, v and w be the velocities 

in the x, y and z directions, respectively, { the downward displacement from equi- 
librium, p the pressure, and p the density of the water. We take 
foro<cz<h 
fork < z< (2) 
=f; forthe 
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Fic. 1.—Longitudinal cross-section of the lake showing the notation. 


Fic. 2.—The vertical equilibrium density distribution in the lake, showing 
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where (; = (2 at z = Ay, and (2 = (3 at z = hy +h, and 


for (z= 0) <2 <m+h(z [top layer) 
= Ug, = Uz, W = We, = pr, p = pr 
forty + =m) < < thetle(z = [middle layer] 


u= ug, U= 3, W= = pa, P = ps 
for hy +he+le(z = m+he)<2z<h [bottom layer] (3) 


where z = {; (z = 0) represents the surface of the water, z = hj +; (z = hy) 
the interface between the top and middle layers, and z = hi +h2+ (z = hi +he) 
the interface between the middle and bottom layers. 

Then, ignoring geostrophic effects, forces of friction, and the viscosity of the 
water, the following equations are satisfied within each layer (Proudman 1953; 
18, 20, 49): 


(i) The equations of motion: 


(ii) The equation of the continuity of volume: 
ox Oz 

(iii) The equation of the continuity of density: 


(6) 


Suffix i takes the values 1, 2 and 3, suffix 1 referring to motion in the top layer, 
suffix 2 to motion in the middie layer, and suffix 3 to motion in the bottom layer. 
The following assumptions are now made: 


(i) Velocities in the y direction may be neglected, i.e. v4 = 0. 
(ii) The vertical acceleration @w;/@t may be neglected in (4). 
(iii) ug, w; and %; are small quantities such that terms involving their squares 


and products may be neglected. Thus in (4) we may neglect the terms 
uy Ox, wy Oz, ug Ow;/ Ox and w; Ow;/ 


| 

Ox oy dz pi Ox 
du, 1 Opi 

Cw; Ow; Ou; Ow; 1 Opi 

(5) 

ops Ops opi 

— + + + =0. 

ot ox ey dz 


Seiches in a narrow lake, uniformly stratified in three layers 
Under these assumptions (4), (5) and (6) reduce to: 


(7) 
(8) 
(9) 


(10) 


In addition, the vertical velocity is given in terms of downward displacement 
by the equation 


(11) 


Equations (7)-(11) are the basic equations to be used in the subsequent analysis. 


5. Solution of the equations 
From (10), (11), 
Opi 


ot Ox ot oz (12) 


Integrating this equation, neglecting products of small quantities, and remember- 
ing that p; (i = 1, 2, 3) must reduce to the equilibrium densities given by (1) when 
ti = 0, yields: 


he a 


a a a3 

To the same order of approximation, it follows from (13) that p; = po at 
z = (z = Ay), and p2 = pg at z = hy le (z = +he), which means 
that density is a continuous function of the depth z during the motion. 


Following the same procedure as Proudman (1953, 361) a solution of equations 
(7){11) is now sought in which 


t 
F("(z)sin sin (14) 


—this denoting a periodic motion starting from rest, in which u; = 0 at x = 0 
and x = /, of wavelength 2//n and period T. F(z) is some function of z, 
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os (15) 


and integrating this gives 
nn, ant 
From (11), (16), 


Integrating this gives 
If (18) is substituted into (13), and then (13) into (8), it follows that 

2ila, l 

hy 

a 


(19) 


Certain boundary conditions on velocity, pressure, and vertical displacement 
have now to be satisfied. If p_ denotes atmospheric pressure, then the boundary 
condition on the pressure at the surface of the water is 


pe at z= % (2 = 0). (20) 


At the horizontal bottom of the lake vertical velocity and displacement must be 
zero, and therefore 


at z=h. (21) 


It is required that vertical displacement be a continuous function of the depth, for 
which 


at z= 
{3 at z= hy +he. 


Velocity and pressure must be continuous across the interfaces between the layers 
so that 


(22) 


= U2, = We, pi= pe at z= =m), (23) 
ug = Us, W2= p2= ps at z= hy = 
To the order of approximation given by assumption (iv) of Section 2, (20) may be 
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written: pi = pa—pogti at z = o, and the continuity conditions on u, w and p, 
given by (23), may be evaluated on the equilibrium interfaces: z = h, and 
2z = hi+hz. Hence, from (20)-(23), we must have: 
(i) pi = pPa-pogti at 
(ii) pi=po at z=hy 
(iii) pe= ps at z= h+he 
(iv) uw at z= hy 
(v) we=ug at hy+he 
(vi) m=m,G =f at 
(vii) we = mg, le = at z= hy+he 
(viii) ws=%3=0 at z=h. 
With these conditions, u, w, ¢ and p derived from the theory are defined for 


2< hy, hy < < +he, and < < h. 
The integration of (19) yields 


art 
fi = pa + =| _ 
(z hy)? 
242 2a, 2lag 
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so that for 24(i), 24(ii), 24(iii) to be satisfied, 


~F\(0) = Fi'(0), 
a 
I I 
—F,(h,) = —F2(h), 
a2 
+he) = + he). 
a2 a3 


From (14), it follows that 24{iv), 24{v) are satisfied if 
= 
+he) = + he), 
while from (16) (18): 24(vi), 24(vii), 24(viii) are satisfied if 
= F2'(hn), 
F2'(hy + he) = + he), 
F3'(h) = 0. 
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Differential equations for Fj, F2 and F3 are now obtained. Substituting p; 
given by (13), into (7), ignoring dt (products of small 
quantities, 7 < 7) and smaller terms, yields 


oui opi 
(34) 


—which, on substituting for u; and p; given by (14) and (25), reduces to 


= (35) 
where 


n*T*g 
og? = 6 
(36) 
The general solution of the differential equation (35) is 
F(z) = Ajcos + Bysin yz (37) 


where A;, B; (i = 1, 2, 3) are constants. Conditions (27)-(33) on F;, when com- 
bined with (37), yield five linear equations in these constants. When the equations 
are solved, expressions for A;, B; are obtained involving an arbitrary multiplying 
constant C,. Substituting these expressions into (37) yields: 


F(z) = aghg cos aghe — x2 Sin aghe sin cos{a;(hy — z)} — 


[ag cos aghg sin + a2 Cos sin aghg] sin{a(; — 2)} 


= cos{xe(h + he—2)}— sin aghg sin{ae(hy + he — 2) 


C. 
Fala) = —cosfas(h—2)}. (38) 
To determine the possible values of «, %2, «3 and period 7, we substitute 
F,(z) given by (38) into (26) and, for convenience put 
ay = (39-1) 
so that from (36), 


ag = gé/he, a3 = ré/hg, (39.2) 
where 


q = = (40) 
The substitution yields the equation 
‘tan(A+£) = (42) 
where 
tan ré + (a2/a3)* tan 


(42) 
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If the positive roots of (41), arranged in ascending order of magnitude, are denoted 
by € = & (k = 1, 2, 3,..., to 00) then, from (39), 


=€/m, = r&k/hs, (43) 
and from (36), 


T = 2m€;/bn (44) 
where 


b = (45) 


To each value of k there are an infinite number of possible modes of oscilla- 
tion of the water having, respectively, wavelengths of 2//n (m = 1, 2, 3,..-, to 0). 
In each mode the period of oscillation T is given by (44), and u, w, ( and p are 
given by (14), (16), (18) and (25), respectively, where F;(z) is given by (38) and 
a1, a2, a3 by (43). Combining these equations yields: 


nbt 
u = sin sin — 
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Cnx is an arbitrary constant and 
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Summarizing the results of the above analysis, possible modes of oscillation 
of the water in the lake have been derived. Equation (41) has an infinite number 
of positive roots which, arranged in ascending order of magnitude, have been 
denoted by £ (k = 1, 2, 3,..., to 00). Corresponding to each root there is an infi- 
nite number of possible modes of oscillation having, respectively, wavelengths of 
2l/n(n = 1, 2, 3,-.-, ©). Thus to every pair of positive integers (k, m) there corres- 
ponds a possible mode of oscillation of wavelength 2//n. Period T, u, w, { and p 
in such a mode are given by (44) and (46). Modes of wavelength 2/ (m = 1) are 
uninodal, those of wavelength 2//2 = 1 (n = 2) binodal, etc. 


6. Ordinary and internal seiches 


When € is small, from (42), tan A ~ (ho+3)&/h; and therefore € tan(A + £) 
~h€*/h;, so that (41) reduces to 


& = hy2/hay. (50) 


Now Ay?/ha; is small, and therefore, from (50), the smallest positive root of (41) is 
given approximately by 


£1 = hy/(hay)'. (51) 


Since A/a; is small the higher roots of (41), i.e. £2, €3,..., are, approximately, the 
roots of 


tan(A+£) = o. 


Consider first the seiches corresponding to £ = £). 
From (44), (51), 


(52) 


T = {2m/bn}{hi/(hai)*} = 2/n(gh). (53) 


Also, ignoring the terms of order £,2, it follows from (48) that Z(z, £1) = (h—2)/(has)*, 
and d{Z(z, £:)}/dz = —1/(hag), so that from (47), 
£1) = — £1)/nm(gh)! = (h—2)/(has)* 


u(z,£1) = £1)/po(gh)* = —(g/as)*. Ge) 
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From (46), (51), (53), sai 


& 2nt 


ant 


o<2<h. (55) 


Formulae (53), (55) may be recognized as those for ordinary seiches. (From 
(55), w and ¢ are greatest at the surface z = o, and the relations for ordinary seiches 
between u, w, p and the surface elevation: —{(z = 0), given by Proudman (1953, 
220-228) are satisfied. Also, (53) reduces to Merian’s formula when m = 1.) 


The smallest positive root £ = £; of (41) gives, therefore, ordinary surface seiches 
in the lake. 


Next consider the seiches corresponding to £ = éo, é3,..., given by (52). 
Eliminating A from (52), using (42), yields 


cos é[cos gf sin rf + (a2/a3)* sin gé cos ré] 
= sin [sin gf sin rf — (a2/ag)* cos g€ cos }(a/a2)*, (56) 
or, alternatively, 
[(a2/a3)* + (a1/ag)! + + 1) sin{(g+r+ 1)£}+ 
+ [(@2/a3)* — (a;/ag)* — (a,/a2)* + 1] sin{(g+r—1)£}+ 
+ [(a2/a3)* + (a1/a3)* — — 1] sin{(g—1+ 1)€}+ 
+ [(@2/a3)* — (a;/a3)* + — 1] sinf(g—r—1)£} = 0. (57) 


Equation (57) may be solved numerically for € (e.g. using a graphical method fol- 
lowed by the Newton-Raphson process as described by Hartree (1955, 190-195)) 
yielding £ = £ (k = 2, 3,...). The seiches corresponding to each value of k have 
periods given by (44). In each seiche mode u, w, ¢ and p are given by (46) where 


u, ®, ¢ and p, derived from (47), (48), (56) may be expressed as follows: 
= —(E&/nb)w(z, Ex) 


= [cos és + (2) sin ré 
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cos cos — (2) sin gf sin 
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con +ho- cos — 


(2) si in sin ré| 
for hy < < 


for <2 <h. (59) 


It is clear from (58) that { = @ = o when z = 0, so that from (46) { = w =o 
when z = 0 for all x and t, i.e. the surface elevation of the water remains zero even 
though the water is displaced vertically below the surface. The seiches corres- 
ponding to the roots = £ (k = 2, 3,...), are, therefore, internal seiches in the 
lake. Those of wavelength 2//n will have periods which are considerably greater 
than the period of the surface seiche of the same wavelength, since £9, §3,..., 
determined from (52) will, in general, be much greater than £ determined from 
(51). 

The above results show that in an internal seiche the surface elevation is zero. 
This is true to the order of approximation in the present theory: Mortimer has 
pointed out during a discussion with the author that in reality it will not be zero 
but will be small in comparison with internal vertical displacements. 


7. Conditions for the validity of the neglect of vertical acceleration 


The preceding theory has been developed after neglecting vertical accelera- 
tion éw/at in the equation of motion, 


(60) 
From (46), the theory gives 
ow 
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Substituting these — into the right-hand side of (60) yields 
nbpoto(z, £x)/Ex 
d{p(z, 
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The neglect of éw/ét in (60) is therefore justified by the results of the theory, 
given by (61), if, from (62), |A| < 1, which, using (47), (48), may be written: 


€1 (@=1,2,3). (64) 


The argument leading to this condition breaks down for ordinary seiches (k = 1) 
since d{p(z, €1)}/dz = o from (54). For ordinary seiches, following the procedure 
given by Proudman (1953, 223), it may be shown that the neglect of vertical 
acceleration in (60) is justified by the results of the theory, insofar as this neglect 
does not significantly affect the pressure , if 


n22h2 
(65) 


The inequalities (64), (65) are necessary conditions for the validity of the neglect 
of vertical acceleration; (64) applies to internal seiches, and (65) to ordinary 
seiches. 


8. Seiche motion in a narrow lake 


Observations of temperature in Lake Windermere (Mortimer 1952a) have 
shown that a gale blowing down the length of the lake displaces the water verti- 
cally, as well as horizontally, as shown by the vertical movements of isotherms 
in a longitudinal section of the lake. When the wind forces maintaining the 
displacement of the water drop to zero longitudinal seiche oscillations originate 
from the displaced condition (assumed to be an equilibrium one). The oscilla- 
tions continue, damped down by friction, until the next wind sets up further 
disturbances. 

During the calm period immediately following the gale, the motion in the lake 
may be regarded as consisting of the sum of the ordinary and internal seiches 
derived in Section 6. Taking all the possible values of k and n these are infinite 
in number. If however ordinary seiches are ignored, together with the internal 
seiches corresponding tom > N andk > K (where N and K are positive integers) 
then the motion in the lake is, from (46), given by 
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The circumstances which govern the choice of N and K are described below. 

A knowledge of the vertical density distribution in the lake (a mean distribu- 
tion derived from temperature observations in the lake may be taken) is sufficient 
to determine u, w, { and p, given by (66), apart from the constants Cy. The 
method of determination is to fit a theoretical density distribution of the type 
given by (1) to the density distribution in the lake (e.g. Figure 4). The parameters 
of this theoretical distribution, together with the results derived in Section 6, 
then determine the roots ¢ and the functions %, #, { and p for substitution 
into (66). 

The modal characteristics of the seiches, considered separately, may be 
examined without having to determine the various values of Cy,x (e.g. Figure 
5 (a)}{h)). However, to determine completely the motion of the water, given by 
(66), these constants have to be found. They may be determined from a knowledge 
of the distribution of temperature in a longitudinal section of the lake when the 
wind forces drop to zero (e.g. Mortimer 1952a, Figure 8b). Isotherms may be 
drawn for this distribution giving a picture of the vertical displacements of water 
in the lake at the commencement (t = 0) of the seiche motion. Starting conditions 
for the motion are thus obtained which are sufficient to determine the various 
values of Cyx as now described. 
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Fic. 3.—Diagram showing the displacement of an isotherm from its 
horizontal equilibrium position, in a longitudinal section of the lake. 


The displacement of an isotherm is shown diagrammatically in Figure 3. 
Letting {ox-(x) denote the initial downward displacement of the isotherm from 
its horizontal position at a depth z = 2%, we have: 


= for z= zy, t=0. 
Substituting this condition into (66) gives 
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If the variation in depth of the isotherm (at t = 0) is z = z0x(x), from Figure 3, 


and therefore 


= sor(x) dx. (69) 
0 


Knowing the temperature distribution in a longitudinal section of the lake at = 0, 
isotherms may be drawn with various mean depths z = z,. Their respective 
variations in depth may be represented by z = z0x-(x), k’ = 1, 2,..., K’, where 
K’ is the number of isotherms required to give a good representation of the 
displacement condition of the water at t = 0. Knowing 2ox-(x), 2% and Cox-(x) 
may be obtained for each isotherm using (68) and (69). Taking k’ = 1, 2,..., K’, 
(67) then gives a set of K’ equations in the constants Cy,. 


To determine Cx from these equations, fox-(x) is represented by a cosine series 
of N terms: 


N 
n=1 


== cos de. (71) 
0 


In (70), N is chosen large enough (while satisfying (64)) for the cosine series to 
give an adequate representation of [ox-(x) for each value of k’. From (67), (70), 


K — 
> = (72) 
k=2 


Since k’ = 1, 2,..., K’, (72) constitutes, for each value of n, a set of K’ linear 
equations in the K—1 unknown constants C,x. Under these conditions, providing 
that K—1 < K’, estimates of the values of C,,% may be obtained using the method 
of least squares (Whittaker & Robinson 1944). 

The condition K—1 < K’ imposes an upper limit on K. Having fixed the 
value of K’ as described above, the best theoretical representation of the seiche 
motion is therefore likely to be obtained by taking K—1 = K’, for this choice 
includes in the analysis the maximum number of different seiches, of nodality n, 
for which the constants C,,% can be determined. 


9. Seiche modes for Lake Windermere 


Mortimer (1952a) made a detailed study of isotherm movements in Lake 
Windermere northern basin from 1947 June 10 to June 13, during a calm spell 
of weather following a north-westerly gale. The course of the internal seiche 
oscillations was followed by taking hourly observations of the vertical distribu- 
tion of temperature at a station situated near one end of the lake. Also, the 
distribution of temperature in a longitudinal section of the lake was observed on 
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June to at the beginning of the observations, and on June 13 at the end of the 
investigation. 

From the observations a mean vertical temperature distribution in the lake for 
the period of observation was derived, from which a curve of mean vertical density 
distribution was deduced (Mortimer 1952a, Figure 15) making the assumption 

that density depended entirely on temperature. Working from this distribution 

: the present theory is now used to calculate the seiche modes (and their periods) 
in the lake corresponding to k = 1, 2,..., 8 (Section 6). 

A density distribution of the type given by (1) is fitted to the experimentally 

determined distribution as shown in Figure 4. From the former, reading off 
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Fic. 4.—Mean vertical density distribution in Windermere northern i os 
basin, 1947 June 10-13. To this is fitted a density distribution of the 


type postulated in the present paper, also shown in the figure. 
p is in g/cm’. 
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values of depth and density (indicated symbolically in Figure 2) we get: 


54m = 0-00012 g/cm$ 
he = poht2/ae = 00007 g/cm 
hg = 31°6m pohs/az = 0-000055 g/cm* 
po = 09991 g/cm 
and from (40): g = 2-3241, r = 1°6377. These numerical values in (57) give 
3°9013 sin(4-9618£) — 1-6783 sin(2-9618¢)— 3-186 sin(1-6864£) — 
— 1°3416sin(o-31365£) = o. 
Solving this equation numerically, the first seven roots = & (k = 2, 3,..., 8) 
have been determined; also = from (51). The results are as follows: 
k I 2 3 4 5 6 7 8 
0700393 0°39293 -1°3535 1°9247 2°5179 3°2412 3°7994 4°2140 
Taking / = 6-6km (for Windermere northern basin) and g = 981 cm/s?, 
periods of oscillation in the seiche modes corresponding to k = 1, 2,..., 8; 
nm = I, 2,...,'7, have been determined from (44) and are presented in Table 1. 


u and { giving the variations with depth of horizontal velocity u and vertical dis- 
placement { in the modes, have also been obtained (from (54), (58), (59)) and are 


Table 1 


Seiche periods (in hours) for Lake Windermere northern basin, 
calculated from the present theory 


n=1 n=2 n=3 n=4 n=5 n=6 


o-181 00905 0036 0°030 
18-1 9°05 3°62 
88-5 44°2 17°7 
149 74°5 29°8 
175 87°5 35°0 


on Aut wn 


plotted out in Figure 5 (a)-(h). The values of u and ¢ are determined, for each 
seiche, to within an arbitrary multiplying constant Cyx (equation (46)). The 
various values of Cy, for the period 1947 June 10-13 have not been determined, 
using the method given in Section 8, since the initial temperature distribution 
corresponding to the state from which the seiche oscillations during the period 
originated (Mortimer 1952a, Figure 8b) is incomplete in the lower depths of the 
lake. It has not been possible therefore to determine vertical displacements of 
the water in the lake for comparison with experimental observations. 
For the surface seiches (k = 1, m = 1, 2,..., 7): 


< = 0-0098, 
and for the internal seiches (k = 2, 3,.... 8; m = I, 2,..., 7): 
< < < = 0°027, 
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i.e. n*n2h?/2I2 and n®n%h;2a;/l2£,2a, do not exceed 0:03. Conditions (64), (65) for 
the validity of the neglect of vertical acceleration are therefore satisfied ignoring 
errors of 3 per cent and less. 

Figure 5(a) gives the variations of ¢ and u with depth z in the ordinary seiches 
(k = 1): at any time instant, { is greatest at the surface of the water and u does not 
vary with depth z. 

Figures 5(b)-(h) give the variations of ¢ and u with depth z in the internal 
seiches (k = 2, 3,..., 8). For the seiches of differing nodality corresponding to 
k = 2 (Figure 5(b)) ¢ is at a maximum at z = 9°5m, i.e. in the lower part of the 
thermocline region. The uninodal seiche, for which T = 18-1 hours from Table 1, 
clearly corresponds to the uninodal internal seiche determined by Mortimer 
(1952a) from a simple two-layered model of the lake. In this model the interface 
between the layers was taken at various depths between 6 and 9-5 m and the uni- 
nodal seiche periods obtained varied between 18 and 19-3h. For the seiches of 
differing nodality corresponding to k = 3 (Figure 5(c)) ¢ has maximum values at 
depths of 5-5 and 20m, its greatest value being at the latter depth. The uninodal 
seiche, for which T = 62-2h from Table 1, would appear to correspond to the 
uninodal internal seiche associated with the lower interface of the three-layered 
model of the lake considered by Mortimer (1952a). Periods calculated for the 
latter seiche varied between 34 and 57h. The seiches corresponding to k = 4 
(Figure 5(d)) and k = 7 (Figure 5(g)) are remarkable for the high values of ¢ 
below the thermocline (i.e. at depths greater than, say, 10m) as compared with 
low values at and above the thermocline. 


10. Conclusions 


From the results obtained by applying the present theory to determine pos- 
sible periods and modes of longitudinal seiche oscillation in Lake Windermere, it 
is concluded that: 


(i) There are internal seiches in which the vertical displacement of the water 
is greatest in the region of the thermocline (Figure 5(b)). These seiches corres- 
pond to those which may be derived from a simple two-layered model of the lake. 
Oscillations of temperature observed at the thermocline appear to be mainly con- 
trolled by these seiches. 


(ii) There are internal seiches of longer period in which the vertical displace- 
ment of the water is greatest below the thermocline in the lower depths of the lake 
(Figures 5(c, d, g, h)); in other seiches vertical displacement in the lower depths 
attains values comparable in magnitude to its greatest values at, and above, the 
thermocline (Figures 5(e, f)). Because of these characteristics, it seems likely 
that the presence of the seiches might be able to explain the large vertical oscilla- 
tions of temperature observed below the thermocline in the lower depths of the 
lake (Mortimer 1952a). Internal waves of this kind are likely to give rise to con- 
siderable horizontal water movements in the lower depths. Mortimer (1951, 
1952a) has suggested that this kind of movement could explain the turbulent 
diffusion of heat and dissolved salts (Mortimer 1941, 1942) in the bottom waters 
of various lakes during thermal stratification. 


(iii) A feature of the internal seiche modes are reversals of horizontal velocity 
as depth increases from the surface to the bottom of the lake in a vertical line. The 
number of reversals increases as the order of the mode increases (i.e. as k increases). 
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Fic. 5. (a)-(d).—ii(z, in and {(z, for Windermere northern 
basin, 1947 June 10-13, giving the variations of u and { with depth z 
in the seiches corresponding to k = 1, 2, 3, 4. At levels denoted by x 
there is a discontinuity in di/dz. 

(a) k = 1, T = 0°181/n hr (ordinary seiches). 

(b) k = 2, T = 18+1/nhr. 

(c) k = 3, T = 62°2/nhr. 

(d) k = 4, T = 88-5/nhr. 
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basin, 1947 June 10-13, giving the variations of u and { with depth z in 
the seiches corresponding to k = 5, 6, 7, 8. At levels denoted by 
there is a discontinuity in du/dz. 

(e) k = 5, T = 116/nhr. 

(f) k = 6, T = 149/n hr. 

(g) k = 7, T = 175/nhr. 

(h) k = 8, T = 194/nhr. 
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For the higher-order modes, the reversals tend to be concentrated in the thermo- 
cline region. Flow reversals of this kind have been observed in experiments with 
model tanks (Mortimer 1952b). 

(iv) To determine the extent to which the various seiche modes influence the 
motion of the water in the lake, particularly in its lower depths, further detailed 
experimental information on water movements below the thermocline is required. 
This information might take the form of either observations of temperature 
oscillations, measurements of water velocities (if this is possible for velocities of 
not more than a few centimetres per second), or ‘‘time—position” records of drift 
bodies in the lake. Mathematically the problem is the determination of the 
constants C;,,; in equation (66). Before actual displacements and velocities of water 
can be calculated the values of these constants must be known. 

(v) The determination of horizontal velocity in the various seiches is of interest 
from the point of view of examining the stability of the internal waves. An 
examination of Figures 5(c)-(h) suggests that vertical velocity gradients in the 
thermocline region might be quite steep, possibly leading to flow instability 
there. Against this, it has to be remembered that vertical density gradient is 
greatest in the thermocline region, and that velocities are likely to be small in the 
seiches of higher order and longer period. 
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On Thermodynamics of Planets* 


C. Lomnitz 
(Received 1960 October 18) 


Summary 

A thermodynamical theory of planets is outlined. General results 
by Prigogine, de Groot and others are applied to the case of a steady- 
state system with a constant temperature distribution. An “excited 
state” is defined by pressure perturbation introduced by an earth- 
quake at the surface of the planet. The resulting transient flows 
towards the perturbed region are analysed, and it is shown that the 
energy transient is logarithmic in time. 

A condition of realizability of the earthquake problem is defined. 
For a given planet type there is a critical size below which no seismic 
activity can occur on the planet. 


Notation 


chemical affinity 

mass concentration of chemical 

phase k 

energy per unit volume 

total energy 

surface heat flow rate per unit 
area 

external force per unit mass on 
phase k 

chemical reaction rate 

material flow of phase k with 
respect to the centre of 
gravity 

heat flow rate per unit area 

Onsager’s phenomenological co- 
efficients 

pressure 

pressure perturbation (due to 
earthquake) 

initial value of AP 


q radioactive heat supplied per 


unit mass 
surface of planet 
specific entropy of phase k 


time 

temperature 

energy per unit mass 

volume per unit mass 

specific volume of phase k 

velocity of centre of mass 

velocity of centre of mass of 
phase k 

total volume of planet 

conjugated force corresponding 
to Jk 

conjugated force corresponding 
to Ja 

coefficients in equations (24) 
and (27) 

chemical potential per unit 
mass of phase k 

density 

density of phase k 

entropy production rate per 

unit volume 

stoichiometric coefficient of 
component k. 
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x. All planets are open thermodynamic systems inasmuch as they may exchange 
energy with the surrounding space. Consider a planet having an indefinite num- 
ber of state variables. In general the following basic equations will apply to the 
system (de Groot 1952a): 


Conservation of mass equation: 
=-p div v. (1) 
Force equation: 
grad P+ > Fxpx. (2) 
t k 
Energy equation: 
2 
= div(Pv+Jo)+ Fr. (3) 
Gibbs’ equation (Second Law of Thermodynamics): 
ds du dv dcx, 


An entropy balance equation may be obtained from the above by elimination 
of the kinetic energy of the centre of gravity between (2) and (3), followed by intro- 
duction of (1) and (4): 


ds 
pT— = — divJg+ -Je+ pare. (5) 
dt k 
This equation can be re-stated as follows, after de Groot: 
ds 
= —divJ,+¢ (6) 


which says that the change of specific entropy in the system is due to the negative 
divergence of an entropy flow J;: 


Js = Ja- 2 (7) 
plus an entropy production o: 
o = (Jog. Xut+ 2 Je. Xe + AJe)/T- (8) 


The vectors X, and X, and the scalar A are called “‘forces” conjugated with 
the flows Jg, Jz and J-. From (5) we deduce that the values of these forces are 


Xu = —(grad 7)/T (9) 
= Fy—Tgrad(ux/T) (10) 
A=- 2 (11) 


2. Up to now a general formulation of the system has been drawn up, closely 
following de Groot. Let us now assume that the temperature T' at every point of 
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the system remains constant in time. The distribution of T in space is unspecified, 
as is the distribution of all other variables of state. We know that if one variable 


is fixed with respect to time a steady state of the first order is obtained, and all 
flows except the flow corresponding to the fixed variable vanish: 


Jk = Je =0 (12) 
a = = (Jg. Xy)/T. (13) 


The latter is the steady-state value of the entropy production, which is also 
its minimum value for that particular temperature distribution. Therefore, if a 
perturbation in any state variable other than temperature is introduced the entropy 
production will increase temporarily. This will be called an “excited state’, to 
be designated by an “x” superscript. Steady state parameters will likewise be 
distinguished by the superscript “‘o’’. 

Let us assume a pressure perturbation AP at every point in the system. Such a 
perturbation might be the result of an earthquake altering the stress distribution 
in a region of the surface of the planet. The excited value of the force X; (equation 
(10)) may then be written: 


x 0 
= X,°- (grad grad). (14) 


The chemical potential 4 is expressed as a function of the state variables as 
follows (de Groot 1952b, p. 113): 


gradp, = grad T+ grad P+ > r,p grad c. (15) 
i 


If all state variables except P remain practically constant we may write (14) thus: 
Xx? = X,°— vz grad AP. (16) 


3- Prigogine (1947) has shown that the excited state is described by the appearance 
of transient flows directed in such a manner as to counter the perturbation and to 
re-establish steady-state conditions. As an example let us consider a transient 
flow of matter J,” defined by a phenomenological equation of Onsager (1931): 


= 2 + (17) 


As we already know, in the steady state Jx° = o since all flows other than the 
heat flow vanish. As the temperature remains unchanged we may write Xy* = X,,° 
and equation (17) becomes 


or, from (16): 
= 2 Lieve grad AP. (19) 


Examination of the entropy equation (8) shows that all coefficients Liz must be 
positive, since otherwise the entropy production o could become negative in contra- 
diction to the Second Law of Thermodynamics. Therefore one can always write 


Ludgrad AP)? > o (20) 
which in turn leads to (cf. equation (19)): 
grad AP > o. (21) 
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The inequality above indicates that the transient flow of matter J,” is always 
of the same sign as the gradient of the perturbation AP. Similar considerations 
apply to electrical, chemical or any other transient flows induced by the perturba- 
tion. 

To sum up, if an earthquake reduces the stress level in a region there will be 
transient flows of matter, electricity, towards the perturbed region. This conclu- 
sion is at least not contradicted by experience; in fact it would be easy to see a 
causal connection between earthquakes and some forms of volcanism, magmatic 
intrusions and orogeny in general. 


4- The development of these transient processes with respect to time is of par- 
ticular interest. Consider the entropy production in the excited state: 


The term in brackets represents the increase in entropy production due to 
the flow of matter, and other transient flows. If we continue to consider flow of 
matter only we may write: 


of = of +[ (23) 


Prigogine has shown that the entropy production in the excited state is always 
of the following general form: 


a2 
o%(t) = + + + (24) 


in which the instant of the perturbation is taken as ¢ = 0. In general, and par- 
ticularly for long times, one may write with sufficient approximation: 


o%(t) = 0° +a/(t+8). (25) 


Comparing this with (23) and designating by AP; the initial value of AP at 
t=0o: 


grad AP = grad AP;[y/(t+8)}* 


1/y = 


(26) 
(27) 


The energy per unit volume expended at any point in the system during a 
time interval ¢ is given by: 


© 


t 
e(t) = dt. (28) 
0 


Let us now introduce (23) and (26) and carry out the integration. The result 
is: 


e(t) = ( 2 AP;)?[log(t + 8) — log 8}. (29) 
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Finally, the total energy expended in reducing the perturbation during the 
time interval t can be found by integrating (29) over the volume V of the planet: 


E(t) = | ( Luvn*)(grad AP; aV. (30) 


If we let t > 00 we realize that E must equal the energy of the earthquake plus 
the total energy of the aftershock sequence. The form of E indicates that the 
aftershock process should be logarithmic in time which is indeed the case (Benioff 
1951). An explanation of aftershocks as caused by local stress adjustments due to 
the influx of energy E, is a natural consequence of this theory. 


5. At the other extreme, one may also view the seismic process as an intermittent 
regulating device, somewhat like a series of steam safety valves on a boiler. In 
this case the temperature can no longer be assumed constant, since the time scale 
is now much ampler than before. 

This perspective is useful chiefly in order to clarify the criteria of realizability 
of the problem. Let g be the heat production per unit volume, f the surface heat 
flow and e the earthquake energy production. Then we have the inequality 


e<qV-fS (31) 


> fs (32) 
or, since V/S = r/3 in a sphere: 


which implies 


q > 3flir, (33) 


this being a necessary but not sufficient condition of realizability. 

If, in the lifetime of a planet, its heat production falls below the amount 
required by (33), all seismic activity on the planet will cease. Conversely, if we 
consider a series of planets of different sizes, all having similar heat productions 
and heat flow rates, there will be a critical radius below which no earthquakes can 
occur. 

The surface of such planets will be relatively smooth and free from the scars 


of orogeny. They will be devoid of atmosphere, oceans and other traces of past or 
present volcanic activity. 


Instituto de Geofisica y Sismologia, 
Universidad de Chile, 
Santiago, 
Chile: 
1960 October. 


References 

Benioff, H., 1951. Bull. Seismol. Soc. Amer., 41, 31. 

de Groot, S. R., 1952a. Thermodynamics of Irreversible Processes. (North- 
Holland Publishing Company (Amsterdam).) See notation at beginning of 
this paper. 

de Groot, S. R., 1952b. Thermodynamics of Irreversible Processes, p. 113. (North- 
Holland Publishing Company (Amsterdam). ) 

Onsager, L., 1931. Phys. Rev., 37, 405. 

Prigogine, I., 1947. Etude thermodynamique des processus irréversibles. (Paris: 
Dunod.) 


; 
. 
ey. 
: 
3 
| 


Comments and Some Suggestions on Hunter’s 
Formula of Reduction of Observed Values of Gravity 
to the Earth Model for use in Stokes’s Integral 


J. C. Bhattacharji 


(Received 1960 October 31)* 


Summary 

The reduction formula as evolved by Hunter, for reducing the 
observed values of gravity to the surface of the Earth model, is discussed 
and proved to be not quite accurate. The main defect involved 
appears to be due to his ignoring altogether the gravity effect on 
account of the sphericity of the Earth while dealing with the difference 
of attraction between the actual and model Earths. 

The corrected reduction formula posing again a delicate problem 
for geodesists, a suggestion is here made for a more convenient method 
of reduction based on a suitable isostatic hypothesis for use in the 
Stokes’s integral without any real difficulty. 


1. Hunter’s method of reduction for Agz 


In Figure 1, let plm represent the actual Earth and PLM the Hunter’s Earth 
model with smoothed topography having slopes not exceeding 0-01, derived from 
weighted averages of heights A of the actual land topography or — 0-615 times the 
depths d of the actual ocean beds over circular areas of about 100 miles radii. 
Again let pab and PAB represent the surfaces of the two hypothetical plateaux, 
the former being identical with the equipotential surface of the actual Earth, that is, 
the geop Gp as defined by Hunter, through p a point of height A» on the actual 
Earth and the latter with the equipotential surface of the Earth model, that is, the 
geop Gp as defined by Hunter, through P the corresponding point of height Hp 
on the Earth model. 

The portion indicated by horizontal shading in the same diagram then repre- 
sents the departure of topography of the actual Earth p/m from the plateau pad 
and the appropriate gravity effect usually known as the terrain correction at p 
due to the actual removal of the above departure of topography is denoted by O; 
when considered all round the Earth and by O when considered up to the distance 
of about 100 miles radius; while that due to the density layer p.i~hp), that is, 
the condensation of the above departure of topography, on the surface of the 
plateau pab is —O,+O when considered all round the Earth and only zero when 
considered up to the distance of about 100 miles radius. Similarly the portion 
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shown by vertical shading in the same figure represents the departure of topo- 
graphy of the Earth model PLM from the plateau PAB and the appropriate 
gravity effect or the terrain correction at P due to the actual removal of the above 
departure of topography is denoted by 9), when considered all round the Earth 
and by Q when considered up to the distance of about 100 miles radius; while 
that due to the density layer p(H~ Hp), that is, the condensation of the above 
departure of topography, on the surface of the plateau PAB is —Q)+Q when 
considered all round the Earth and only zero when considered up to the distance 
of about 100 miles radius. 


Let gp be the observed value of gravity at a point p of height A» on the surface 
of the actual Earth plm, g,' the corresponding value of gravity at p after the actual 
removal of the departure of topography of the actual Earth plm from the plateau 
pab and g,’’ the corresponding value of gravity at p after condensation of the above 
departure of topography, on the plateau pab so that 


= &pt+O1 (1) 
and 


= &p = gp t+ O1—O1+0 = gpt+ O. (2) 


Similarly let gz be the reduced value of gravity at the corresponding point P 
of height Hp on the surface of the Earth model, gg’ the corresponding value of 
gravity at P after the actual removal of the departure of topography of the Earth 
model PLM from the plateau PAB and gz" the corresponding value of gravity at P 
after condensation of the above departure of topography, on the plateau PAB 
so that 


ge = get+Qy (3) 


ge” = gr’ Q = = get Q. (4) 


Now let us introduce a density layer — ps(h— Ay) on the plateau pab and assume 
the surface of the plateau pab to be equivalent to an infinite plane, tangent to the 
surface of the plateau at p. Then the gravity effect at p due to the above density 
layer on the infinite plane surface becomes zero instead of —(—O,+0O) when 
the sphericity of the Earth’s plateau is taken into consideration. 
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Let us introduce another density layer +p.H—/y) forming a doublet with 
the former and equivalent to two separate density layers +p(H—Hp) and 
—pdhp—Hp) on the plateau PAB and assume as before the surface of the 
plateau PAB to be equivalent to an infinite plane, tangent to the surface of the 
plateau at P. Then the gravity effects due to the above density layers on the 
infinite plane surface become separately zero instead of —Q)+ 2 and — A(hy— Hp) 
respectively, A being the factor for the usual mass reduction, and that due to the 
removal of the layer —p<(hpy—Hp) from the infinite plane surface becomes 
—A(hy—Hp) instead of —2A(hp—Hp) when the sphericity of the Earth’s 
plateau is taken into consideration in both the cases. 

Thus the value of gravity gz’’ at P on the surface of the Earth model with 


its departure of topography from the plateau PAB condensed on the surface of 
the plateau reduces to: 


when the sphericity of the Earth’s plateau is neglected 


&p" = O+ F(hp—Hp)— A(hy— Hp) (5) 
or when the sphericity of the Earth’s plateau is considered 
ge" = Flhp— Hp)—2A(hyp— Hp)—- (6) 


where F denotes the factor for the usual free-air reduction. 


Consequently the value of gravity gz at P on the surface of the Earth model 
with its topography in place becomes: 


when the sphericity of the Earth’s plateau is neglected 


8B = O+F(hy—Hp)— A(hp—Hp)-2 (7) 
or when the sphericity of the Earth’s plateau is considered 
ge = £pt+O1+ (8) 


Now let yp be the theoretical value of gravity at P on the associated spheroidal 
surface of the same potential as that at P on the Earth model. Then the gravity 
anomaly Agg as derived by Hunter at P on the surface of the Earth model by 
ignoring the effect of the sphericity of the Earth’s plateaux pab and PAB becomes: 

Age = &p+O-Q— A(hp—Hp)+ F(hp—Hp)—yp 
= Hp)]+ A(hy — Hp) 
= A(hy — Hp) 
= Ag+O-Q- A(hy— Hp) (9) 


where Ag denotes the usual free-air gravity anomaly at p on the surface of the | 


actual Earth. 


But the assumption that the gravity effect on account of the sphericity of the 
Earth or its plateaux is practically negligible does not seem to be quite accurate. 
For as already indicated in relations (6) and (8) when the sphericity of the Earth’s 
plateaux pab and PAB are considered, the gravity anomaly Ag, at P on the surface 
of the Earth model with its topography in place, becomes 


Age = gp+ O1—Q)—2A(hp—Hp)+ F(hp—Hp)—yp 
= + O1 —2A(hp—Hp) 
= 


(10) 
| 
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The formula (10) is obviously very different from the formula (g) derived by 
Hunter. 


2. An alternative method of reduction for Agz 


The reduction formula (10) can be derived more easily as follows: 

Let gp be as before the observed value of gravity at a point p of height Ay on 
the surface of the actual Earth plm and gy’ the corresponding value of gravity at p 
after the actual removal of the departure of topography of the actual Earth p/m 
from the plateau pab so that 


&p = &p+O1. 

Let gz’ be the value of gravity at the corresponding point P of height Hp 
on the surface of the plateau PAB, reduced from the value of gravity gp’ at p on 
the surface of the plateau pab by actually removing or adding the shell pP accord- 
ing as hy is greater or less than Hp and gg the corresponding value of gravity at P 


on the Earth model after restoring the smoothed topography of the Earth model 
in place so that 


ge’ = + F(hp—Hp) —2A(hp—Hp) 
= &p+Oi+ F(hp—Hp)—2A(hp — Hp) 
ge = ge 
= gpt+O,+ 
Now let yp be as before the theoretical value of gravity at P on the associated 
spheroidal surface of the same potential as that at P on the Earth model. Then 
the gravity anomaly Agg at P on the surface of the Earth model with its smoothed 


topography in place and without neglecting the effect of the sphericity of the Earth’s 
plateaux pab and PAB, reduces to 


Age = gp+O1—Q)+ F(hp —Hp)—2A(hp—Hp)—yp 
= &p—[yp— F(hp—Hp)]+ O1— 21 — 2A (hy — Hp) 
= 8p—ypt O1— —2A(hp— Hp) 
= Hp). 


Thus the value of Agg as derived by Hunter differs from that as obtained now 
by the amount: 


— A(hp — Hp) —O+Q. (11) 


It now remains to be seen whether the effect: O,;— Q,— A(hp—Hp)-O+Q 
is really significant at all points on the Earth’s actual surface. 

From the very definition of the Earth model it appears obvious that whatever 
be the nature of the Earth’s actual topography, the average heights of the indivi- 
dual compartments of gravity zones at least beyond 300 miles radius, are not likely 
to differ appreciably from those of the Earth model while the gravity effect of the 
Earth’s actual topography is comparatively small. Keeping the above in view let 
us examine the significance of the quantity: O,—Q);—A(hp—Hp)—O+Q in 
the following three different cases under ideal conditions: 
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Case (i) 
We have from Figure 1 on substituting the values of O, and Q; 
— A(hp — Hp)—- 
= [0+ A(1+x)(hp —Hp)]—[Q+ Ax(hyp —Hp)|-—O+ Q— A(hy—Hp) 
since the effects of the portions of the shells of thickness (1+.x)(hp—H>) lying 
at a distance of about 100 miles radius, are respectively 
A(1+x)(hp~Hp) and Ax(hp—Hp) = o. (12) 


This shows that so long as the average heights of zones Hayford 18 to 1 for the 
actual Earth, having, except at most for only a few nearer zones, practically the 
same value as that for the Earth model, remain at or below the level of the equipo- 
tential suface PAB, the relations (10) and (9) are practically identical. 


Case (ii) 
We have from Figure 2 on substituting the values of O; and Q; 
— A(hp— Hp) 
= [0+ Ax(hp— Hp)]—[Q+ A(1 —x)(hp Hp)]|—O+ Q- A(hy— Hp), 


since the effects of the portions of the shells of thickness x(hp—Hp) and 
(t—x)(Ap—Hp) lying at a distance of about 100 miles radius are respectively 


Ax(hp—Hp) and A(1 —x)(hp— Hp) 


= —2A(1—x)(hp—Hp) (13) 
which does not vanish unless x is equal to 1. 


Thus it appears that so long as the average heights of a number of compart- 
ments of zones 18 to 1 for the actual Earth, having, except at most for only a few 
nearer zones, practically the same value as that for the Earth model, remain above 
the level of the equipotential surface PAB of the Earth model, the expression (13) 
is likely to be an appreciable quantity and becomes a maximum when x reduces 


to zero. 
Case (iti) 
We have from Figure 3 on substituting the values of O; and Q) 
0, — A(hp— Hp) 
= [0+ Ax(hy—Hp)]—[Q+ A(1+x)(hp— Hp)]—O+ A(hy— Hp) 
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for reasons similar to case (i) 


= —2A(hp—Hp) (14) 


which does not vanish unless hp and Hp become identical. 

This shows that so long as the average heights of a number of compartments 
of zones 18 to 1 for the actual Earth, having, except at most for only a few nearer 
zones, practically the same value as that for the Earth model, remain at or above 
the level of the equipotential surface pab of the actual Earth, the expression (14) 
is likely to be an appreciable quantity and remains constant for all values of x. 

It is to be noted in this connection that when p happens to be below P, we 
again get the similar results provided proper care is taken in the case of signs. 
Thus the Hunter’s reduction formula (9) seems to be in error by an amount 
equal to O, — Q; — A(hy— Hp)—O+Q which, as is evident from the above, does 
not reduce to a negligible quantity except when the Earth’s actual topography is 
likely to fulfil the conditions enumerated in case (i) considered above. Though the 
given cases are only imaginary, yet the results are likely to give a clear indication of 
the order of error that may occur because of the approximations involved in the 
reduction formula (9) derived by Hunter. 


Actual Earth 
Earth Model 


3. Practical difficulties of Hunter’s model Earth anomalies 


Since an accurate evaluation of the quantity O,;—Q, is never an easy matter 
even with the best available gravity tables, because of the uncertainties involved in 
the estimated values of average heights of distant zones, obtained from maps 
however accurate, the reduction formula (10) as derived by the author, though quite 
correct from theoretical considerations, proves to be really difficult in its practical 
application. Again from relation (10) it becomes clear that Agg may quite correctly 
be regarded as a type of Bouguer anomaly reduced to the surface through P of the 
Earth model instead of the geoid, after having considered the topography out to the 
antipodes instead of the specified distance of only about 100 miles from the point. 
But due to the fact that Agg is a Bouguer anomaly modified as above, its reduction 
does involve mass-transfers within and beyond the limit of about 100 miles radius, 
which ultimately affect the equipotential through P of the actual Earth, giving 
rise to a new equipotential surface or better called a co-geop corresponding to the 
Earth model and of the same potential as that of the equipotential surface PAB 
of the Earth model. Moreover the resulting vertical separation between the above 
equipotential surface through P of the actual Earth and the corresponding co-geop 
due only to the topographic effect of the mass-transfers, though normally quite 
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small except in regions of rugged topography, is, as already discussed in my pre- 
vious note (1959), quite difficult (not exactly impossible as stated in my previous 
note) to be accurately obtained even after heavy reductions. It is, however, true 
that since the average heights of zones beyond about 400 miles radius for the actual 
Earth, are practically equal to those for the Earth model, the topography beyond 
about 400 miles radius does not appreciably affect the potential at P. But since the 
estimated average heights of zones 18 to 12 are normally likely to be in error by 
150m or so particularly in the more distant ones, the order of uncertainty involved 
in the computed values of the vertical separations between the equipotential surface 
through P of the actual Earth and the co-geop, on that account alone, can easily 
be up to 2m which is never negligible. This also is definitely a disadvantage in 
the accurate computation of the deviations of the equipotential surface of the 
actual Earth through P. Thus even though one is prepared to take all possible 
care to carry out rigorous computations of O; — Q, in order to reduce the gravity 
anomaly to the surface of the Earth model by making use of the relation (10) and 
also make proper allowance for the deviations of the equipotential suface of the 
actual Earth through P due to the mass-transfers involved in the reduction, it 
would still be very difficult to reduce them with the order of accuracy actually 
needed in order to apply them correctly in the Stokes’s integral. 


4- Suggestion for a new type of Earth model anomalies Agz< 


In the circumstances, to secure the advantages of an isostatic reduction in order 
not only to diminish the topographic effect of the mass-transfers on both the 
gravity and the potential but also to obtain the result with more ease and precision, 
“the obvious way is to make the reduction on a frankly isostatic basis, choose the 
most appropriate depth of compensation and frankly face the fact that the transfer 
of matter as we do in imagination in such a reduction, means also a warping”’ of 
the equipotential surface through P of the actual Earth, though mostly insignificant 
being considerably smaller in magnitude than what is generally obtained in the 
absence of the effect of any isostatic compensation. Even if the effect on the devia- 
tion of the equipotential surface becomes sometimes more significant, allowance 
may easily be made accordingly to reduce the gravity anomaly to the level of the 
terroid (according to Dr Hunter’s notation (de Graaff-Hunter 1958)) of the same 
separation as the equipotential surface through P of the actual Earth has from the 
co-geop, which is essential in problems connected with the figure of the Earth. 
Moreover the corresponding shift of the centre of mass and of the rotation axis 
due to the mass-transfers implied in isostatic reductions, is also likely to be ex- 
tremely negligible, though it does not affect the final isostatic anomaly reduced to 
the level of the terroid of the same separation as the equipotential surface through P 
of the actual Earth has from the co-geop on taking the indirect effect of gravity into 
consideration as stated above. Having thus known the gravity field on the surface 
of the terroid, we can easily determine the elevations and depressions of the co-geop 
from the associated spheriodal surface by the application of the Stoke’s integral. 

But obviously these are not the final results sought. What are wanted are the 
elevations and depressions of the equipotential surface through P of the actual 
Earth from the associated spheroidal surface. Lambert & Darling (1936) have 
designed necessary tables for computations of the indirect effects of gravity as well 
as the corresponding deviations of the equipotential surface through P of the actual 
Earth. Now there being also no sensible shifting of the centre of mass and of the 
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rotation axis, the computed values of the deviation of the equipotential surfaces 
through P of the actual Earth may be straightaway used in deriving the elevations 
and depressions of the equipotential surface through P of the actual Earth from the 
associated spheroidal surface without any further considerations. It may be em- 
phasized here that it is quite immaterial for the present problem to consider 
whether the adopted isostatic hypothesis holds good for the actual Earth or not; for 
in any case we do obtain an elevation or depression of a co-geop above the corres- 
ponding equipotential surface which requires the gravity anomaly to be actually 
reduced to the level of the terroid in order to be correctly applied to the Stokes’s 
integral giving the elevation or depression of the co-geop from the associated 
spheroidal surface which again when added to the corresponding vertical separa- 
tion between the co-geop and the equipotential surface through P of the actual 
Earth gives the required elevation or depression of the equipotential surface 
through P of the actual Earth from the associated spheroidal surface. In other 
words the co-geops act, in the present problem, more or less as intermediaries to 
take us simply from ground levels to terroids, then from co-geops to associated 
spheroidal surfaces and finally from associated spheroidal surfaces to the equipo- 
tential surfaces of the actual Earth. 

Hence we need not apprehend any real setback on account of our making use 
of the isostatic hypothesis. 

Let us then reduce the gravity anomaly: Ag, isostatically (preferably Pratt- 
Hayford hypothesis) to the surface of the terroid of the same separation 5N as 
the equipotential surface through P of the actual Earth has from the correspondin g 
co-geop and then denote it by Aggc so that 


Agec = &p+O-A(hp— Hp) + He +8N)—yep—Q+AEc+ AEric¢ 
= gp—[yep—F(hp—Hp)+ A(hp—Hp)—O]+F 
= Agat+ (15) 


where AE¢ stands for the effect of compensation for the differences of average 
heights of zones A to O for the actual Earth and the Earth model, which is prac- 
tically the same as the effect of compensation for the heights (h—h’) of zones A 
to O, h being the average height of zones A to O for the actual Earth and h’ the 
average heights of zones A to O for the Earth model; AEz,¢ is the combined 
effect of topography and compensation for the differences of average heights of 
zones 18 to 13 for the actual Earth and the Earth model and is comparatively 
small; the corresponding differences of average heights of zones beyond 13 being 
too small, the combined effect of zones 12 to 1 is negligible; F.5N stands for the 
indirect effect of gravity. 

The relation (15) is obviously much simpler and more accurately determinable 
than the relation (10). The quantities O; and Q); have advantageously been re- 
placed by the easily computable ones O and Q, the latter being generally negligible. 
The quantities Ac and AE7,c being approximately linear functions of the differ- 
ences of average heights of zones for the actual Earth and the Earth model are also 
easily and accurately computable. Last of all, 5N being computed on the hypothe- 
sis of complete isostatic compensation using the same differences of average heights 
of zones A to O only, the corresponding combined effect of zones beyond O being 
negligible, is not only small but also extremely easy to compute or simply read 
from a chart, easily designed for the purpose. 
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Apart from the facilities in actual computations as enumerated above, the 
Earth model isostatic anomalies have other advantages too. Though the anomalies 
are actually based on the isostatic hypothesis, the mass-transfers involved in the 
reductions are notably less than those postulated in the current isostatic reductions. 
Moreover these refer to points on the surface of the terroid closely resembling the 
smoothed surface of the Earth model with comparatively small but accurately 
known vertical separation 5N between them and may thus be termed as near- 
surface anomalies also. Being again freed of the effects of irregular topography 
these are not only more regular and representative of their neighbourhood than 
the usual free-air anomalies, but also comparable with similar anomalies anywhere 
in the world without recourse to any significant assumptions of internal densities 
of the Earth, and as such are more suitably applied in the Stokes’s integral, 
providing ultimately a means of obtaining more precisely the required vertical 
separations at the near-surface points only instead of the geoid, between the 
associated spheroid level and the surface points of the actual Earth. 

Hence the isostatic gravity anomaly Aggc reduced as above, is not only free 
from the main defects of the purely topographic reductions and the large-scale 
mass-transfers involved in the current isostatic reductions, but also quite regular 
and representative, and applies correctly in the Stokes’s integral for the purpose 
of determining ultimately the shape of the actual Earth instead of the geoid, 
without any real difficulty. 
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New Gravity Measurements in West Pakistan 


K. Helbig and H. I. S. Thirlaway 
(Received 1960 November 9)* 


Summary 


The thermostat fitted to the Worden “Master” Gravity Meter 
results in much improved drift characteristics, and gravity differences 
under difficult conditions can now be measured with a precision of 
o-2mgal. This paper describes its use in West Pakistan to establish a 
base network of 8 air links connected to Woollard’s Karachi airport 
station where g is taken to be 978-963 0cm/s?. A selection of 21 new 
values based on these air links is reported, and the stations described. 
They are regularly distributed over the Northern half of the province. 
The gravity air links relative to the adopted value of Karachi are 
estimated to have a standard deviation of 0-2 mgal or better. Standard 
deviations, based on repeat observations by road and rail transport, 
reach a maximum of 0-3 mgal. 

Comparisons between 12 Survey of India pendulum stations and 
gravity meter values are graphed. The Cambridge pendulum obser- 
vations lie within + 2 mgal of the new gravity meter results over a range 
of 7oomgal. It may be possible to derive a sway correction for the pre- 
Cambridge pendulum values from the results. 


1. Introduction 


A Worden “Master’”’ Gravity Meter (No. 551) has been used to establish new 
gravity bases, and to compare old ones, in West Pakistan. The novel feature of 
this model is the provision of a thermostat to assist existing temperature com- 
pensation devices. 

The base network was established at Quetta, Multan, Lahore, Sarghoda, Rawal- 
pindi and Peshawar using the internal services of Pakistan International Airways. 
Quetta and Lahore were linked directly with the National Gravity Base at Karachi 
Civil Airport. Quetta—Lahore—Karachi—Quetta were measured on a round trip 
of 34 hours by Thirlaway, the Lahore—Karachi difference being determined by 
extrapolating the Karachi overnight drift curve. Lahore—Rawalpindi-Sarghoda— 
Lahore—Multan—Lahore were measured in that order within 10 hours. Mr Aziz 
of the Department of Geology, Punjab University, measured the Lahore-Multan 
difference, and Helbig the remainder of the network. 

Quetta and Lahore have also been connected separately, by road and rail, by 
Thirlaway. Intermediate stations in the Punjab, North West Frontier Province, 
and Swat State were measured on these occasions and repeated on the return rail 
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journey. An attempt by Helbig to make a second road link was partially frustrated 
when transport broke down at Mughal Kot for 5 days, and temperature control 
was lost. This event, however, served to emphasize the remarkable consistency 
of drift while the thermostat was in operation. 


Figure 1 locates the stations where gravity differences were observed in this 
series. 
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Fic. 1.—Location map of new stations. 


2. Drift and meter constants 


Table 1 summarizes the rail, road and air connections between Quetta and 
Lahore, and serves to demonstrate the new standards in regularity of drift which 
have been achieved with the use of a thermostat on a Worden Meter. 


Table 1 


Comparison of three independent gravity differences 
Quetta (Geophysical Institute)—Lahore (Department of Geology) 


Date of measurement, 1960 


Mean drift 
1st Quetta 2nd Quetta Lahore 


Transport (mgal/day) 


12 February 6February—-_ Rail 
11 February 


Ag (meal) 


533°2+0°75 
(Mean of 6 readings 
in 6 days at Lahore) 


5 February 


{1523 hours) 


15 March 


13 April 
{0827 hours) 


31 March 26 March- Road 


28 March 


14 April 13 April Air 


(1850 hours) (1145 hours) 


$33'2+0-01 
(Mean of 3 readings 
in 3 days at Lahore) 


533°2 
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In all three connections, drift was determined from the two Quetta measure- 
ments. A tare, generated by manipulation of the large dial, was observed at 
Lahore after the rail journey. Poor large dial technique, and possibly immature 
springs, was responsible. If the large dial is repeatedly brought to its mark clock- 
wise through one revolution, an agreement of 0-1 mgal can be obtained for succes- 


sive readings within a minute or two. Other independent observations so far 
made with this meter are summarized in Table 2. 


Table 2 
Independent gravity values observed with Worden “Master” 
gravity meter No. 551 


g Karachi = 978-963 0 cm/s® 


Date Transport Observer Value of gravity 
(1960) (cm/s*) 


16 March Road Thirlaway 979°3305 


25 May Road Helbig 3302 
16 March Road Thirlaway +1639 


26 May Road Helbig 163 6 


7 May Air “4680 
Sarghoda Helbig 


Mughal Kot 


23 May Road 4679 

The first Tank and Mughal Kot values were measured during the Quetta-— 
Lahore road connection, when the closing drift check was made after an interval 
of 16 days (Table 1). The second measurements at these stations, and the second 
Sarghoda value, are based on the average drift curve (1-62 mgal/day from 44 
observations) for the 17 days ending 1960 May 22 at Lahore. This circumstance 
arose out of the transport failure which caused loss of temperature control on 
the May 27 at Mughal Kot. The second Mughal Kot value, therefore, depends 
on a linear drift extrapolation of 4 days from an observed curve of 17 days’ dura- 
tion. In our experience, the mean drift of this meter is constant, both in transport 
and whilst stationary, until the temperature setting is altered, or temperature 
control is temporarily lost. The scatter of points during the 17 days observations 
was largely due to the tidal effect, and an eye fit was easily made with good precision. 

The meter constants are briefly summarized for future reference. From 1960 
January to 1st March, the thermostat was in continuous operation at go °F. 
Subsequently it was operated at 110°F. These figures are correct to +1 °F. 
The small dial range is 200mgal with a constant of 0-2547(5) mgal/division, 
while the mean large dial constant varies between 8-970 mgal/division at 100 to 
8-988 mgal/division at 800. Above 300 the mean variation is linear. All values 
reported have been calculated with a large dial constant of 8-984 mgal/division. 


3. Results at new stations 


The air links are connected to Woollard’s Station at the Civil Airport, Karachi. 
Values of g determined at this station are given in Table 3. 
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Table 3 
Gravity meter values at Karachi Civil Airport 
Observer 
Year (and reference) Meter g (cm/s?) 
1948 Woollard (1952) Worden 1ob 978-963 4 
1950 Muckenfuss (1952) Worden roc -962 6 
ts Stahl (1958) Western 42EQ 963 7 
1954 Marussi (Stahl 1958) Worden 6 9592 


A nominal value of 978-963 ocm/s? for Karachi has been selected from these 
results. 

Table 4 gives a distributed selection of new stations at which gravity differences 
from Karachi have been measured via the air links. Details of the air connections 
with Karachi are also tabulated. With the exceptions of Chiniot and Rabwah, the 
rail and road links were made in 1960 February and March respectively (‘Table 1). 
Chiniot and Rabwah were measured during the attempted 1960 May road con- 
nection Lahore—Quetta, which failed to close at Quetta due to loss of temperature 
control. These two values are, therefore, based on an extrapolation by 1 day of 
the previous 17 days drift at Lahore (cf. Table 2, the second Sarghoda value 
measured on the same day, and the second Mughal Kot value 3 days later). 


Table 4 


Values at new Pakistan gravity stations based on gravity meter connections 
by air with Karachi Civil Airport, where g = 978-963 0 cm/s? 


Station Description Value of gravity 


cm/s? 
Air Links 
Quetta Geophysical Institute Geomagnetic Block 978°857 5 
Quetta Airport P.1.A. Traffic Office 978-880 7 
Lahore Geology Department Geophysical “tomb” 979°390 7 
Lahore Airport (Walton) Concrete block, verandah in front 979°392 0 
P.I.A. Office 
Sarghoda Airport Front of VIP room entrance 979°468 o 
Rawalpindi Airport (Chaklala) Left end verandah facing P.I.A. Office 979°3500 
Peshawar Airport Left end verandah facing P.I.A. Office 979°395 0 
Multan Airport Verandah opposite booking office door 979°255 3 
Rail Links 
Railway Stations, Platform level, at Booking Office 
Latitude Longitude Height 
° , ” ° , ” ft 
Mach 29 51 55 67 19 30 3246 978-762 8 
Jacobabad 28 16 35 68 26 oo 183 979°191 4 
Rahim Yar Khan 28 25 oo 7o 18 05 979°143 7 
Khanpur 28 37 36 70 39 00 295 979°133 4 
Bahawalpur 29 24 05 71 39 'I5 380 979°200 4 
Khanewal 30 18 35 71 55 50 430 979°283 3 


Montgomery 39 39 10 73 06 550 979°3247 
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Table 4—continued 
Road Links 
Latitude Longitude 

opposite MES Inspection 30 67 978880 4 

bungalow 
Manikawa. Fort Gate 31 69 
Mughal Kot. Hospital entrance 31 70 
Bannu. Mission bungalow 

Tochi Road, entrance 33 7° 
Peshawar. 38 Dean’s Hotel, 

entrance 34 71 979°397 4 
Saidu Sharif. 8 Swat Hotel, 

entrance 34 72 979°305 9 
Mansehra. Dak bungalow, 

entrance 34 73 979°245 2 
Murree. Survey of Pakistan, 

entrance 33 73 
Mangla. Dam Rest House, 

entrance 33 73 979°371 8 
Dera Ghazi Khan. PWD 

Inspec. Bungalow entrance 3° 47 70 979°197 6 
Rawalpindi. Christchurch, 

Primary B.M. 
Multan Cantt. Church, Identical with Survey of Pakistan 

Primary B.M. New Series (1960) 
Chiniot. PWD Insp. bungalow 

B.M. 
Rabwah. T.I. College B.M. 


979°099 5 
979°163 9 


979°347 


979°0440 


979°349 6 
979°247 0 


979°4617 
979°4733 


4- Comparisons with other measurements 


Whenever possible, Survey of India pendulum stations were re-occupied 
during the course of new measurements (Survey of India 1926, 1931 and 1936). 
The pendulum results are compared with those from the gravity meter in Table 5 


and in Figure 2. Karachi and Rawalpindi gravity meter values are corrected for 
large location differences. 


Table 5 


Comparison of Survey of India pendulum results, with gravity meter 
results relative to Karachi = 978-963 0 cm/s® 


g (cm/s*) 


Pendulum station 


Pendulum Gravity meter 

Qila Saifullah 978°970 -969 7 
Fort Sandeman 979°062 063 3 
Tank 979°329 "3305 
Nowshera 979°421 "4207 
Sibi 979°120 "1196 
Quetta 978°851 ‘8537 
Multan 979°243 "246 
Montgomery 979°321 "3253 
Bahawalpur 979°199 "199 8 
Khanpur 979°132 "1345 
Karachi 978961 959 

Rawalpindi 979°346 "3471 
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Fic. 2.—Comparison of pendulum and gravity meter measurements 

in West Pakistan. The five points represented by a combined cross 

and solid circle compare pendulum observations with Worden meter 
No. 102. 
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Figure 2 also illustrates comparisons with the Survey of Pakistan Worden 
gravity meter No. 102. Pendulum stations in Kalat and Mekran were observed 
with this meter by Thirlaway. In this paper they have been corrected by 
—2-omgal from published values (Thirlaway 1960) because the new Quetta- 
Karachi air link difference has reduced the value formerly used at the Quetta 
base by this figure. The remaining results from Meter No. 102 plotted in Figure 
2 were observed by Dr G. Norgaard in 1953 using the Karachi base. 


5- Discussion 

The air links have established a maximum difference of 587-2 mgal between 
West Pakistan airports served by Pakistan International Airlines. These dif- 
ferences are useful for checking changes in gravity meter calibration constants 
used in Pakistan. All measurements were completed with a temperature con- 
trolled Worden meter, “Master’’ model, No. 551 within 6 months of its calibra- 
tion by the manufacturers using their standard geodetic meter practice on pendu- 
lum bases near Houston. 

The base network of air links has been used to obtain gravity values at new 
stations in N. Baluchistan, the Frontier Province and Punjab areas of West 
Pakistan, and to compare values at pendulum stations and other gravity meter 
stations. The pendulum values are not inconsistent with the value of 978-9630 
cm/s? adopted from gravity meter values established at Karachi Civil Airport. 
Over the small pendulum range of 700 mgal it is not possible to identify a calibra- 
tion error in Meter No. 551. 

Repeat measurements at stations measured with Survey of Pakistan Worden 
Meter No. 102 show a systematic error due almost certainly to wrong calibration 
constants. On the advice of the manufacturers, Meter No. 551 is to be assumed 
correct and an error of 0-4 per cent presumed in the calibration factor of Meter 102. 
Independent repeat observations under difficult circumstances with Meter No. 551 
suggest that, using the thermostat, an internal precision of 0-3 mgal is possible 
on large dial work, with drift checks at intervals of several days. Use of the 
small dial only, under similar conditions, would give a precision of 0-2 mgal 
while daily drift checks bring a precision of 0-1 mgal within easy compass on the 
small dial. This new meter, in conjunction with a base network of air links, 
removes the need for tedious and expensive loops. For non-believers in the use 
of the large dial, considerable savings can be effected by specifying a small dial 
range of 200 mgal or even 400 mgal in place of the standard 80 mgal. 
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Reports of Discussions 


Geophysical Discussions of the Royal Astronomical 
Society 


A Geophysical Discussion entitled The Geological Time Scale, was held on 
1961 January 27, under the chairmanship of Professor L. R. Wager, who in his 
introduction noted that the discussion was to centre upon the validity of the data: 
that is to consider the experimental errors and the natural perturbing factors in 
radioactive age determinations. After describing the initial suggestion of the helium 
method by Rutherford and the pioneer contributions of Rayleigh, Holmes and 
Boltwood, Professor Wager drew attention to the critical improvement in instru- 
mentation, which the advent of Neir’s mass spectrometer represented. Thus 
modern radioactive age determinations virtually begin in 1937. After a period 
in which little work was done in this country the groups at Harwell, Cambridge 
and Oxford have in the last decade made major contributions and it is from these 
groups that the speakers come. 

The Chairman then introduced the first speaker Mr J. E. T Horne, who des- 
cribed the lead method of dating uranium and thorium minerals. He regretted 
that there was not sufficient time to consider the interesting, but as yet imprecise, 
dating of common lead minerals such as galena. The uranium and thorium 
minerals are somewhat restricted to intrusive rocks so that they are of limited value 
for calibrating the geological time scale, but provide reliable dates against which 
to check the other methods. When both uranium and thorium occur in a mineral, 
chemical and lead isotope analyses enable the age to be determined from three 
independent ratios: U238/Pb206, U235/Pb20? and Th232/Pb298, Since the two uranium 
isotopes are present in the same proportions a fourth age though not a further 
independent one, may be deduced from the Pb?°?/Pb26 ratio. To obtain reliable 
ages from these ratios two conditions must be fulfilled. First, there should be no 
lead in the mineral when it formed. If any had been present then it would be 
revealed by the nonradiogenic isotope Pb, which provides a means of correction. 
Second, the mineral must remain a closed chemical system throughout its history. 
Some consequences of the loss or gain of uranium, thorium or their decay products 
were considered by the speaker. Thus it was noted that the Pb?°7/Pb2° ratio is 
only affected if there is differential leaching of either of these isotopes. Again 
loss of radon will not affect the U235/Pb207 ratio. It was stressed that each speci- 
men should be analysed chemically as well as isotopically in order to check the 
various ratios against each other, that more than one specimen should be studied 
from each locality, and that polished or thin sections should be examined for 
signs of weathering. The aim should always be to obtain concordant ages rather 
than to deduce the true age from discordant results. 

The second speaker was Mr J. A. Miller, who described the use he had made 
of the potassium-argon dating method. He first explained that the method 
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depended on the decay of K*° to A*® by K electron capture and that for any such 
system to be suitable for dating purposes, the decay constants must be known, 
none of the daughter product must be initially present and that the daughter 
product once formed must not be lost or augmented. In the K4°/A4° method the 
first two conditions are fulfilled and it is in the loss of the daughter product that 
difficulties arise. Thus errors in the age determination are due to experimental 
errors and to loss of argon. The experimental errors arise in the measurement of 
the isotope ratio of the sample, the measurement of the argon volume and the 
determination of the potassium content. The error in the isotope ratio is small 
compared with the other two. A one per cent error in either of these will give 
approximately a one per cent error in the age. Losses of argon from minerals are 
less readily assessed and not all potassium bearing minerals are suitable. The 
micas have proved the most successful having a high potassium content and high 
argon retentivity. Hornblende and sanidine have also been used successfully. 
Mr Miller next considered the effects of regional metamorphism. Recrystalliza- 
tion may be expected to release some or all of the argon from the lattice. Such a 
process has a definite activation energy, above which the loss would rapidly be 
complete and below which little happens. Results obtained from the Moine 
schists show this effect. These rocks have been metamorphosed more than once 
but the marked consistency of the results suggests that all traces of previous meta- 
morphic episodes have been obliterated and the last recorded. Consequently ages 
obtained from the schists indicate the age of the last metamorphism and shed no 
light on the age of the sediments themselves neither do they give information about 
the earlier metamorphisms. These results agree very closely with those obtained 
from similar material by Dr Moorbath using the Rb/Sr method. The problem 
of thermal metamorphism can be more complex. Experiments have shown that 
argon is lost below 600 °C and that if the temperature is raised to 800°C most 
but not all of it will be removed. Should the biotite now be removed from the 
crucible it will appear perfectly fresh and the X-ray pattern is unchanged from 
the original. Thus it is impossible to tell whether a sample of biotite has lost 
argon by mild thermal metamorphism. The speaker said that there was evidence 
that the argon was held in the mica lattice in two sites, from one of which it was 
easily removed while from the other escape was more difficult. As yet the nature 
of the sites is not clear nor how much argon is assigned to each. Possibly in mild 
thermal metamorphism all the argon from the low energy site may be lost or 
possibly some may remain. Results from near a granite with no distinct thermal 
aureole have given an anomalous age in the country rock 100 yards from the 
contact. More detailed work is being carried out to see if there is a continuous or 
discontinuous transition from the age of the country rock to the age of the intrusion. 
In all the determinations weathered rocks have been avoided because of the ob- 
vious possibility of argon loss. In concluding the speaker commented that using 
the K/A method any changes which take place in the parent mineral result in the 
loss of argon which means that problems of migration are obviated. 

In opening the discussion on this paper the Chairman asked what the importance 
of the time factor was in the loss of argon, to which the speaker replied that once 
the activation energy had been exceeded a slow loss rate was maintained. In reply 
to Mr Dodgson, the rate of loss of argon was compared with normal diffusion 
losses, showing that the large loss at 600°C and 1 200°C separated by low losses 
at intermediate temperatures suggested a more complex process. A discussion of 
the accuracy of the potassium determination with flame photometers then 
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followed, in which Dr Vincent said that for determining potassium at the relatively 
high concentrations encountered in micas and potash felspars, a simple measur- 
ing instrument, using a coal gas~compressed air flame, was probably to be pre- 
ferred to more complex photometers working with higher temperature flames, 
since relatively few elements other than Na and K are excited and much potential 
interference is thus eliminated at source. A series of potassium determinations 
on a given sample should all agree within a standard deviation of about o°5—1 per 
cent of the amount of K2O present. The procedure is so rapid that at least 4 to 
6 replicate determinations can be made on each sample, so increasing the con- 
fidence in the average value obtained. Having performed it very many times, 
Dr Vincent felt that the classical gravimetric method of J. Lawrence Smith was 
both less precise and less accurate, besides being most time-consuming and 
requiring considerable skill. It is essentially an extraction method, and tends to 
yield negative systematic errors, especially with materials so difficult to de- 
compose as the micas. Should the methods for extracting argon become 
sufficiently refined to be applied directly to rocks low in potassium, such as 
basalts, then neither gravimetric nor flame photometric methods were likely to 
be sufficiently precise or accurate. Neutron activation analysis might perhaps 
prove the best method for such materials. Dr N. Snelling said that in some 
rocks, which he had dated, the potassium in the micas had been replaced by 
magnesium to give chlorite, so that the percentage of potassium was reduced 
from 8 to 1 per cent and yet the age determinations were compatible with the 
regional pattern. One therefore wants to know whether these changes take place 
during the initial cooling or gradually later. Mr Miller replied that he had had no 
experience with chloritization, but that when he had examined weathered micas 
it was iron which appeared to be moving. 

After a brief historical introduction to the rubidium-strontium method 
Dr Moorbath showed that the ages of common igneous and metamorphic minerals, 
such as the micas and potassium felspars can be measured with a precision of 
+2-3 per cent by the mass-spectrometric isotope dilution technique. The 
Rb/Sr ages of igneous rocks, which have not suffered disturbing physico-chemical 
effects since their formation, agree closely with the ages obtained by other methods 
and date the time of crystallization. The ages of the Shap and Dartmoor granites 
were discussed as examples of this type. In metamorphic rocks, particularly those 
formed from originally igneous material, different minerals tend to give different 
ages; The pattern which appears is that the muscovite age is greater than the 
potassium felspar age, which is itself greater than the biotite age. The muscovite 
age approaches the age of crystallization, while the biotite age reflects the age of 
the latest metamorphism. Consequently the age patterns can yield information 
of a succession of events which have affected a rock. Several examples from 
recent research were quoted, including measurements from the Lewisian and 
Moinian metamorphic complexes of the Scottish Highlands. Thus the ages of 
the Laxfordian and Scourrie metamorphisms were recorded by the biotites and 
felspars respectively in certain pegmatites from the North West Highlands. In 
contrast, in comparatively fine grained high grade metamorphic rocks, for example 
mica schists of sedimentary origin, the different minerals give concordant Rb/Sr 
ages which are in good agreement with the K/A ages and date the last meta- 
morphism. Since igneous rocks such as granites tend to remain closed systems 
with respect to the overall Rb/Sr ratio during subsequent metamorphism, a whole 
rock Rb/Sr measurement can also yield the age of crystallization of the granite 
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In the discussion which followed this paper Dr T. F. Gaskell said that it was 
most convincing that two different schools at Oxford and Cambridge get the same 
date for the Moine series to within the stated error of measurement. It gives us 
confidence that we are gradually building up a geological time scale, with accu- 
rately fixed markers. The Chairman commented that although the radioactive 
determinations could not do better than stratigraphical palaeontology at its best, 
it was rapidly becoming a real aid to geologists. Dr Snelling added that radioactive 
methods are of supreme importance in precambrian successions in which strati- 
graphical palaeontology is powerless and emphasized that the precambrian repre- 
sented a far longer period of time than the subsequent 5 or 6 hundred My. Dr R. 
St]. Lambert commented on the Geological Time Scale exhibited, saying that the 
accuracy of the subdivisions varied greatly. Divisions of the Tertiary were 
generally well-established in some favourable localities, such as the Western U.S.A., 
where interbedded volcanics contained biotite and sanidine, so that accurate 
measurements were possible. The Cretaceous and Jurassic periods were quite 
well defined, but the Permo-Triassic formations were not dated in any but broad 
outline. The general outlines of the Carboniferous and Devonian were known 
from stratigraphically well-controlled plutonic intrusions, but the Lower Palaeo- 
zoic boundaries were dated on the scantiest of evidence. In particular, there were 
no certain Silurian dates, while the lower parts of the Cambrian were not dated 
at all. Mr Dodgson drew attention to an age of 650 My for a Russian lower Cam- 
brian sylvite by the K/Ca method and Mr Keen commented on the Rb/Sr results 
from M.L.T. on authigenic Cambrian glauconite. Miss C. M. Botley mentioned 
the discordance of the two general determinations of the age of the Earth, the 
one by U2%5 giving 3350 My and the other by the Rb/Sr method 4500 My. She 
added that the age cannot be much greater than the U25 value, because if it were 
U235 would not be found in nature. The Chairman said that the Rb/Sr method 
clearly involved some problems that were hard to understand and that there 
appeared to be a tendency for migrations from larger crystals into the smaller 
accessory minerals such as epidote and apatite. The speaker agreed with this, 
adding that the safest technique was to complete determinations with all the 
minerals likely to carry Rb and with the total rock, and then finally to analyse 
those minerals which do not normally carry Rb. Continuing, he repeated that 
the granites did appear to form closed systems which allowed an important 
application of the age determinations; it is possible to study the distribution of 
dated granites. In North America an approximately concentric pattern has 
been found with the oldest intrusions in the centre of the shield area. More 
recently Gastill has summarized data from other regions and found no similar 
correlation between ages and distribution in terms of present-day continental 
masses; in fact even the North American pattern is complicated by new re- 
sults. Thus the continents have existed at their present margins for at least 
3,000 My, and rejuvenation does not follow any simple pattern, but occurs with 
the restoration of the geosynclinal state by factors which remain largely unknown 
and unpredictable. Mr Miller added that in his Japanese work he had found the 
younger rocks on the continental side of the island arc, which is inconsistent 
with the idea of the outward growth of continents. Dr N. Thorley, from King’s 
College, Newcastle, asked if there were any chance of dating rocks of palaeo- 
magnetic interest. Mr Miller said that some work had been done in conjunction 
with the Rock Magnetism group at Cambridge and other work had been planned 
with Dr Clegg’s group. In reply to Mr Nelson, it was stated that the rhenium/ 
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osmium method had been used in some of the Canadian work but it suffered 
from the extremely limited occurrence of rhenium. 

In closing the discussion Professor Wager said that the development of age 
determination methods is still at a very interesting stage. We are getting more 
reproduceable results and we should, of course, like to improve this. However, 
results by diffetent methods are now known to be not necessarily the same, nor 
do the Rb/Sr ages on whole rocks and separated minerals always give the same 
answers. The correlation between different kinds and intensities of metamorphism 
suffered by the rocks and minerals, which are apparently responsible for the 
difference in ages by the rubidium-strontium and potassium—argon methods, is 
only just beginning to be investigated. It seems that age patterns sometimes record 
complicated series of events in the history of a rock—ranging from mild regional 
or thermal metamorphism to almost complete refusion. 


A Geophysical Discussion was held on 1961 February 24 on the subject of 
Upper Atmosphere Ionization and Aurorae. Mr J. A. Ratcliffe took the chair. 

Dr J. W. Dungey (A.W.R.E.) dealt with the origin of aurora and associated 
phenomena under moderately quiet conditions. Most theoretical study to date 
has followed Chapman and Ferraro in considering the flow of plasma past a mag- 
netic dipole, no magnetic field being assumed in interplanetary space, but Hoyle 
and Alfvén have both discussed the effect of such a field, and it now appears that 
the problem deserves a lot more study. Pioneer V found a steady interplanetary 
field of a few gamma, probably perpendicular to the ecliptic. Dr Dungey argued 
for an approximately southward interplanetary field using ground-based magnetic 
observations. 

A qualitative picture of the flow round the Earth can be obtained from hydro- 
magnetics. There are neutral points near the equatorial plane, when the field is 
southward. The limiting lines of force from the neutral points between them 
cover two surfaces, one of which is a doughnut intersecting the Earth in closed 
curves surrounding the poles, to be identified as instantaneous auroral zones. 
An important feature of this model is that the plasma flow inside the doughnut 
has a general direction opposite to that of the interplanetary plasma. Now the 
electric field generated by the flow is fed along lines of force into the ionosphere, 
and the reversal of flow in the doughnut corresponds to a sudden reversal of the 
ionospheric electric field at the auroral zones. The electric field drives currents 
which produce magnetic variations at the ground; the currents flow nearly parallel 
to the equipotentials because of the Hall effect. The well known reversal at the 
auroral zone provides a striking argument for the postulated southward inter- 
planetary field. The same electric field is also supposed to cause auroral motions. 

The proposed model provides an accelerating mechanism for the auroral 
primaries which can gain energies of many keV near the neutral points. One 
would expect them to arrive near the auroral zones, because they spiral quite 
tightly around lines of force, but to penetrate down to auroral levels their velocities 
must be nearly parallel to the local magnetic field. This facilitates the investigation 
of their motion, but more work is needed. 

Dr S. Evans (Scott Polar Research Institute) spoke about his further work 
relating the movement of visible auroral features to the local magnetic disturbance. 
At Halley Bay the mean movement of individual auroral features for a given hour 
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of the night was found to be related to the deviation of the horizontal component 
of the local magnetic field from the quiet condition for the same hour. The re- 
lationship, in amplitude and direction, implied the effective transport of negative 
charge by the auroral features. Similar conclusions have been reached by other 
workers on the basis of mean diurnal behaviour of radar and visual aurora, and of 
the magnetic elements. Dr Evans has recently analysed auroral movements from 
all-sky photographs taken at Tromso and it is apparent that the association with 
the magnetic disturbance, noted earlier, may arise by chance. The most noticeable 
feature at all stations in both hemispheres, is the average movement towards the 
west in the early part of the night and towards the east after local midnight; the 
period of change-over may cover several hours. The largest component of the hori- 
zontal magnetic disturbance, for auroral zone stations, is the change from positive 
(northerly) deviations in the early part of the night to larger negative deviations 
later. The results from Tromso show that the magnetic variation reversal occurs 
on the average one to two hours before the reversal of auroral movement. Most 
geophysical quantities show some diurnal variation and a direct relationship must 
not be suggested when the only real one is the common position of the sun. 

Thus an attempt has been made to find the relationship for a number of indi- 
vidual hours. When the horizontal component of the magnetic disturbance is 
compared to the auroral movement vector, it is found that the perpendicular com- 
ponent is, on the average, only a little greater than the parallel component and that 
it is not always of the same sign. Therefore, either the regions of auroral luminosity 
cannot be instantaneously associated with nett electric charge, or the movement 
vector used is an inadequate measure of the bulk auroral motion. 

K. D. Cole has considered an auroral arc, assumed to be a slab of highly con- 
ducting material orientated east-west, when acted upon by a neutral wind having 
a north-south component. The vertical component of the magnetic field causes 
an east-west emf to drive irregularities of electron density (and presumably 
luminosity) along the arc. The ratio of the normal movement of the arc to the 
drift velocity of the electrons was found to be equal to the ratio of the Pedersen 
and Hall conductivities, which is always positive and peaks at about 20 in the 
dynamo region. An examination has therefore been made of isolated arcs, having 
bright features which enable the drift velocity to be measured. It is found that the 
normal and the drift velocities are hot approximately proportional to one another 
and further, the sign is frequently opposite to that predicted. The result shows 
that other, and larger, emf’s dominate auroral movements, or the failure of the 
theory to take account of where the current circuits may close has modified the 
result beyond usefulness. 

Dr C. D. Watkins (R.R.S.) reviewed the shapes and movements of auroral 
radar echo regions. Unwin has shown that diffuse echoes are obtained from a thin 
layer, less than 5km thick, at a height of about r10km extending continuously 
polewards for several hundred kilometres. This type of echo persists for hours at a 
time. It seems likely that the layer of ionization is as extensive in an east-west 
direction, though magnetic field aspect effects complicate direct measurements. 
Discrete echoes have individual range extensions of less than 40km and durations 
of minutes, but often exist in groups spread over several hundred kilometres. 
These echoes occur at about the same height as the diffuse echoes but have height 
extensions of 20km or more. 

Observations at U.H.F. have been carried out largely in North America. 
Presnell & others have shown that diffuse echoes are only detected during daylight 
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and arise from an extensive layer at an average height of 110km whereas discrete 
echoes are obtained from a thin band of ionization arranged along a parallel of 
magnetic latitude and extending in height for several tens of kilometres. At 
200 Mc/s, the heights of both types of echo varied by 40km above and below the 
mean height quoted above, with the height spread decreasing as the observing 
frequency is increased. 

Weak auroral ionization has been observed to be arranged in comparatively 
thin arcs extending for some s50okm along a magnetic parallel. Auroral radio 
echoes are generally accompanied by visible displays but no obvious correlations 
between echoes and particular forms has yet been found. 

The movements associated with auroral echoes have been observed in two 
ways, by gross motions of the echo regions, and by doppler measurements of the 
reflected signals. Both methods show the existence of velocities of up to 
3000m/s. The motions of echo regions vary systematically throughout the day, 
with a westward movement before magnetic midnight and eastward thereafter. 
The direction and magnitude of the movements are related to the sign and size 
of the accompanying magnetic disturbances. Velocities deduced by doppler 
methods have values comparable with those measured from the gross motions 
The high values of velocity observed by both methods indicate that the movements 
are not those of the neutral wind. 

Mr W. R. Piggott (R.R.S.) spoke about the relation between abnormal ioniza- 
tion in the lower ionosphere and auroral activity. Day to day changes in both 
phenomena are closely associated, but ionospheric phenomena which most often 
accompany a general increase in auroral activity are usually more frequent at the 
diurnal minimum of auroral occurrence. This is the central problem to be ex- 
plained in the relations between auroral and ionospheric perturbations. The same 
problem arises when magnetic and ionospheric phenomena are compared. 

The most common phenomena observed are the presence of auroral Es or 
retardation Es. The next most common phenomenon is blackout, which indicates 
the presence of considerable ionization in the D region. When the period of maxi- 
mum zenith auroral activity overlaps with the period when Es types ‘a’ and ‘r’ 
are common, the aurora is usually accompanied by Es type ‘a’ or ‘r’. Similarly 
when it overlaps the period of high blackout incidence it is usually coincident 
with a period of blackout. Statistically, the blackout and visible aurora reach 
appreciably lower latitudes than the Es phenomena, though at different times of day. 

The Es phenomena and the D region phenomena, each show a fairly regular 
and consistent statistical variation with latitude, longitude, and time. In both 
cases the abnormal ionization is mainly confined to a limited part of the day and 
a limited zone of the Earth. For any hour, the zone of activity can be defined 
reasonably accurately for blackout phenomena and for the Es phenomena, and the 
arcs of maximum activity open out in a clockwise direction for the latter. This 
suggests that blackout is associated with the incidence of predominantly positive 
particle streams and Es with negative particles. 

There appears to be little doubt that maximum auroral activity at the lower 
latitudes occurs at times roughly half way between the times of maximum Es and 
maximum blackout incidence. At higher latitudes, where these times are more 
widely separated, there are usually two maxima which correspond reasonably with 
the ionospheric evidence. 

The total incidence of Es is much greater than that of blackout and, statistically, 
the most common diurnal changes in auroral incidence are closely allied to those 
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of Es. Thus both phenomena are usually first seen at high latitudes and later 
move towards the equator. 

However, there appears to be a missing factor which modifies the relations 
between visible aurora and the ionospheric phenomena, particularly at the lower 
latitudes. A possible cause might be the presence of local currents at certain hours 
which could increase the brightness of the aurora and make it easier to detect. 

Opening the discussion, Professor V. C. A. Ferraro pointed out that a mechan- 
ism for the acceleration of corpuscular radiation in the vicinity of the Earth is re- 
quired because the Sun-Earth delay times imply a corpuscular velocity almost an 
order of magnitude lower than that required for penetration down to 100km 
altitude in the atmosphere, or that observed in doppler shifted H-« 
from aurorae. This raised the question of what magnitude of interplanetary 
magnetic field is required for Dungey’s acceleration mechanism to operate. 

Dr J. W. Dungey replied that Pioneer V measured two components of a total 
field of 3gamma under quiet conditions, rising to 15 gamma under disturbed 
conditions. It is frequently supposed that the ring current may be identified with 
the outer Van Allen belt, but the flux density has been observed to weaken during 
the main phase of a magnetic storm, whereas the ring current is required to 
strengthen. The neutral points would require to be outside the Van Allen belt, 
both to explain the latitude of the auroral zone, and to ensure that the Van Allen 
particles are trapped in the Earth’s field. A great enhancement of the southward 
interplanetary field alone might be sufficient to account for the main phase. 
Clarifying a further point raised by the Chairman, Dr J]. W. Dungey said that 
positive particles are assumed to be responsible for evening aurora and negative 
particles for the morning, because of the sign of the north-south electric field 
which must act across the auroral zone in order to drive the disturbance currents. 
However, the occurrence of H-«, indicating the presence of incoming protons, 
is not always restricted to the first part of the night. Mr Piggott’s evidence was 
based on the latitude versus local time distribution of blackout and Es, and the 
interpretation of these in terms of Stermer theory. 

Dr L. Thomas said that from spectroscopic evidence Omholt believes that the 
primary radiation is almost entirely electrons; direct rocket observations in an 
aurora support this view and intensity changes in visible aurora have been corre- 
lated with flux changes in the Van Allen belt directly overhead. Professor Ferraro 
pointed out that there will be separation of the particles as they enter the atmo- 
sphere in any case. It should not be assumed that we shall only observe 
the primaries. 

Mr J. A. Ratcliffe pointed out that the radar echoes were not produced by 
Thomson scattering from free electrons, but by scattering from gradients in elec- 
tron density and thus we should expect doppler and range-change observations 
to agree. However, the direction in which anisotropic irregularities in electron 
density move is not the same as that of the individual electrons, nor the same as 
that of the current. A serious problem in general dynamo theory is the choice of 
zero in calculating the disturbance field. Dr Evans avoided this by comparing the 
change in the magnetic vector end-points with the end-points of the movement 
vectors. Mr J. McDowall asked if Dr Evans subtracted the quiet day variation 
from the total disturbance of H in plotting his diagrams; at Halley Bay it is small 
compared to a moderate storm field, but at Tromso it is not. 

In closing, Dr J. W. Dungey asked for IGY records from as many stations as 
possible to be used to derive the instantaneous disturbance current pattern to 
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form a general picture, on a world wide scale. Dr S. Evans said that Dr Dungey’s 
problem had been tackled with pencil and paper with the 1932/3 Polar Year data, 
but much fuller use of computer capabilities could now be made. Instead of 
assuming an overhead line current to account for the horizontal disturbance at 
each available station, a function having (probably) four unknown parameters 
describing current density and direction at each point on a spherical surface 
surrounding the Earth could be inserted. Each ground station would enable three 
unknowns to be removed over a limited area. Considerations of continuity, and 
perhaps symmetry would have to be used to complete the picture but it would 
enable full use to be made of all three components of the disturbance field. Dr 
C. D. Watkins noted that simple relationships between single radar echoes and 
magnetic bays, and similar visual relationships, should not be expected when 
individual correlations between radar echoes and visual aurorae have not been 
found. Mr Piggott remarked that it was unfortunate that there was not time to 
deal with magnetic and ionospheric disturbances following nuclear explosions in 
the upper atmosphere where the source of disturbance was known. In analysing 
IGY data from different stations it should be remembered that the ionosphere 
behaves differently at each place: account should be taken of the known local 
properties and the relation to the geomagnetic poles. 


A Geophysical Discussion on the Geodetic Uses of Artificial Satellites 
was held on 1961 March 24, under the chairmanship of Dr A. H. Cook. In open- 


ing the meeting the Chairman said that the discussion was concerned with the 
ways in which geometrical relations between points on the surface of the Earth 
could be derived from the observation of satellites. The discussion would not 
consider the information which could be obtained about the external gravity field 
of the Earth from the motions of satellites; a subject which had been discussed 
on a number of occasions in the Society. The fundamental geodetic use of satellites 
is their simultaneous observation from stations whose relative positions we wish 
to determine. In this way one may get information about the shape of the surface 
of the Earth and of its external equipotential surfaces, which could not be obtained 
from ordinary geodetic methods confined as they were to land areas. So far no 
satellites have been used for these geodetic purposes. 

The first speaker, Brigadier G. Bomford, outlined the objectives of the work 
and estimated the accuracy with which the observations must be made. The 
major problem is the connection of the reference spheroid across the oceans. 
Thus for example the form of the American and European arbitrary spheroids is 
well known, but they are not on the same reference datum. Radar links already 
exist across the Atlantic, but these only link latitude and longitude leaving the 
distances to the centre of the spheroids indeterminate. In order to reduce two 
surveyed continents to the same reference datum, simultaneous observations 
are made of a satellite from stations on the two continents and hence the 
geometrical relations between the observing stations obtained. The main 
alternative to this method is the indirect method of observing gravity which until 
very recently has in any case been impracticable in surface ships. We know that 
the geoid departs from any spheroid by + 50m, possibly + 100 m, but not by more. 
Since it is these small scale features which we wish to observe, accuracy to within 
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of Es. Thus both phenomena are usually first seen at high latitudes and later 
move towards the equator. 

However, there appears to be a missing factor which modifies the relations 
between visible aurora and the ionospheric phenomena, particularly at the lower 
latitudes. A possible cause might be the presence of local currents at certain hours 
which could increase the brightness of the aurora and make it easier to detect. 

Opening the discussion, Professor V. C. A. Ferraro pointed out that a mechan- 
ism for the acceleration of corpuscular radiation in the vicinity of the Earth is re- 
quired because the Sun—Earth delay times imply a corpuscular velocity almost an 
order of magnitude lower than that required for penetration down to 100okm 
altitude in the atmosphere, or that observed in doppler shifted H-« 
from aurorae. This raised the question of what magnitude of interplanetary 
magnetic field is required for Dungey’s acceleration mechanism to operate. 

Dr J. W. Dungey replied that Pioneer V measured two components of a total 
field of 3gamma under quiet conditions, rising to 15 gamma under disturbed 
conditions. It is frequently supposed that the ring current may be identified with 
the outer Van Allen belt, but the flux density has been observed to weaken during 
the main phase of a magnetic storm, whereas the ring current is required to 
strengthen. The neutral points would require to be outside the Van Allen belt, 
both to explain the latitude of the auroral zone, and to ensure that the Van Allen 
particles are trapped in the Earth’s field. A great enhancement of the southward 
interplanetary field alone might be sufficient to account for the main phase. 
Clarifying a further point raised by the Chairman, Dr J. W. Dungey said that 
positive particles are assumed to be responsible for evening aurora and negative 
particles for the morning, because of the sign of the north-south electric field 
which must act across the auroral zone in order to drive the disturbance currents. 
However, the occurrence of H-«, indicating the presence of incoming protons, 
is not always restricted to the first part of the night. Mr Piggott’s evidence was 
based on the latitude versus local time distribution of blackout and Es, and the 
interpretation of these in terms of Stormer theory. 

Dr L. Thomas said that from spectroscopic evidence Omholt believes that the 
primary radiation is almost entirely electrons; direct rocket observations in an 
aurora support this view and intensity changes in visible aurora have been corre- 
lated with flux changes in the Van Allen belt directly overhead. Professor Ferraro 
pointed out that there will be separation of the particles as they enter the atmo- 
sphere in any case. It should not be assumed that we shall only observe 
the primaries. 

Mr J. A. Ratcliffe pointed out that the radar echoes were not produced by 
Thomson scattering from free electrons, but by scattering from gradients in elec- 
tron density and thus we should expect doppler and range-change observations 
to agree. However, the direction in which anisotropic irregularities in electron 
density move is not the same as that of the individual electrons, nor the same as 
that of the current. A serious problem in general dynamo theory is the choice of 
zero in calculating the disturbance field. Dr Evans avoided this by comparing the 
change in the magnetic vector end-points with the end-points of the movement 
vectors. Mr J. McDowall asked if Dr Evans subtracted the quiet day variation 
from the total disturbance of H in plotting his diagrams; at Halley Bay it is small 
compared to a moderate storm field, but at Tromse it is not. 

In closing, Dr J. W. Dungey asked for IGY records from as many stations as 
possible to be used to derive the instantaneous disturbance current pattern to 
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form a general picture, on a world wide scale. Dr S. Evans said that Dr Dungey’s 
problem had been tackled with pencil and paper with the 1932/3 Polar Year data, 
but much fuller use of computer capabilities could now be made. Instead of 
assuming an overhead line current to account for the horizontal disturbance at 
each available station, a function having (probably) four unknown parameters 
describing current density and direction at each point on a spherical surface 
surrounding the Earth could be inserted. Each ground station would enable three 
unknowns to be removed over a limited area. Considerations of continuity, and 
perhaps symmetry would have to be used to complete the picture but it would 
enable full use to be made of all three components of the disturbance field. Dr 
C. D. Watkins noted that simple relationships between single radar echoes and 
magnetic bays, and similar visual relationships, should not be expected when 
individual correlations between radar echoes and visual aurorae have not been 
found. Mr Piggott remarked that it was unfortunate that there was not time to 
deal with magnetic and ionospheric disturbances following nuclear explosions in 
the upper atmosphere where the source of disturbance was known. In analysing 
IGY data from different stations it should be remembered that the ionosphere 
behaves differently at each place: account should be taken of the known local 
properties and the relation to the geomagnetic poles. 


A Geophysical Discussion on the Geodetic Uses of Artificial Satellites 
was held on 1961 March 24, under the chairmanship of Dr A. H. Cook. In open- 


ing the meeting the Chairman said that the discussion was concerned with the 
ways in which geometrical relations between points on the surface of the Earth 
could be derived from the observation of satellites. The discussion would not 
consider the information which could be obtained about the external gravity field 
of the Earth from the motions of satellites; a subject which had been discussed 
on a number of occasions in the Society. The fundamental geodetic use of satellites 
is their simultaneous observation from stations whose relative positions we wish 
to determine. In this way one may get information about the shape of the surface 
of the Earth and of its external equipotential surfaces, which could not be obtained 
from ordinary geodetic methods confined as they were to land areas. So far no 
satellites have been used for these geodetic purposes. 

The first speaker, Brigadier G. Bomford, outlined the objectives of the work 
and estimated the accuracy with which the observations must be made. The 
major problem is the connection of the reference spheroid across the oceans. 
Thus for example the form of the American and European arbitrary spheroids is 
well known, but they are not on the same reference datum. Radar links already 
exist across the Atlantic, but these only link latitude and longitude leaving the 
distances to the centre of the spheroids indeterminate. In order to reduce two 
surveyed continents to the same reference datum, simultaneous observations 
are made of a satellite from stations on the two continents and hence the 
geometrical relations between the observing stations obtained. The main 
alternative to this method is the indirect method of observing gravity which until 
very recently has in any case been impracticable in surface ships. We know that 
the geoid departs from any spheroid by + 50m, possibly + 100 m, but not by more. 
Since it is these small scale features which we wish to observe, accuracy to within 
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of Es. Thus both phenomena are usually first seen at high latitudes and later 
move towards the equator. 

However, there appears to be a missing factor which modifies the relations 
between visible aurora and the ionospheric phenomena, particularly at the lower 
latitudes. A possible cause might be the presence of local currents at certain hours 
which could increase the brightness of the aurora and make it easier to detect. 

Opening the discussion, Professor V. C. A. Ferraro pointed out that a mechan- 
ism for the acceleration of corpuscular radiation in the vicinity of the Earth is re- 
quired because the Sun-Earth delay times imply a corpuscular velocity almost an 
order of magnitude lower than that required for penetration down to 100km 
altitude in the atmosphere, or that observed in doppler shifted H-« 
from aurorae. This raised the question of what magnitude of interplanetary 
magnetic field is required fur Dungey’s acceleration mechanism to operate. 

Dr J. W. Dungey replied that Pioneer V measured two components of a total 
field of 3gamma under quiet conditions, rising to 15 gamma under disturbed 
conditions. It is frequently supposed that the ring current may be identified with 
the outer Van Allen belt, but the flux density has been observed to weaken during 
the main phase of a magnetic storm, whereas the ring current is required to 
strengthen. The neutral points would require to be outside the Van Allen belt, 
both to explain the latitude of the auroral zone, and to ensure that the Van Allen 
particles are trapped in the Earth’s field. A great enhancement of the southward 
interplanetary field alone might be sufficient to account for the main phase. 
Clarifying a further point raised by the Chairman, Dr J. W. Dungey said that 
positive particles are assumed to be responsible for evening aurora and negative 
particles for the morning, because of the sign of the north-south electric field 
which must act across the auroral zone in order to drive the disturbance currents. 
However, the occurrence of H-«, indicating the presence of incoming protons, 
is not always restricted to the first part of the night. Mr Piggott’s evidence was 
based on the latitude versus local time distribution of blackout and Es, and the 
interpretation of these in terms of Stormer theory. 

Dr L. Thomas said that from spectroscopic evidence Omholt believes that the 
primary radiation is almost entirely electrons; direct rocket observations in an 
aurora support this view and intensity changes in visible aurora have been corre- 
lated with flux changes in the Van Allen belt directly overhead. Professor Ferraro 
pointed out that there will be separation of the particles as they enter the atmo- 
sphere in any case. It should not be assumed that we shall only observe 
the primaries. 

Mr J. A. Ratcliffe pointed out that the radar echoes were not produced by 
Thomson scattering from free electrons, but by scattering from gradients in elec- 
tron density and thus we should expect doppler and range-change observations 
to agree. However, the direction in which anisotropic irregularities in electron 
density move is not the same as that of the individual electrons, nor the same as 
that of the current. A serious problem in general dynamo theory is the choice of 
zero in calculating the disturbance field. Dr Evans avoided this by comparing the 
change in the magnetic vector end-points with the end-points of the movement 
vectors. Mr J. McDowall asked if Dr Evans subtracted the quiet day variation 
from the total disturbance of H in plotting his diagrams; at Halley Bay it is small 
compared to a moderate storm field, but at Tromso it is not. 

In closing, Dr J. W. Dungey asked for IGY records from as many stations as 
possible to be used to derive the instantaneous disturbance current pattern to 
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form a general picture, on a world wide scale. Dr S. Evans said that Dr Dungey’s 
problem had been tackled with pencil and paper with the 1932/3 Polar Year data, 
but much fuller use of computer capabilities could now be made. Instead of 
assuming an overhead line current to account for the horizontal disturbance at 
each available station, a function having (probably) four unknown parameters 
describing current density and direction at each point on a spherical surface 
surrounding the Earth could be inserted. Each ground station would enable three 
unknowns to be removed over a limited area. Considerations of continuity, and 
perhaps symmetry would have to be used to complete the picture but it would 
enable full use to be made of all three components of the disturbance field. Dr 
C. D. Watkins noted that simple relationships between single radar echoes and 
magnetic bays, and similar visual relationships, should not be expected when 
individual correlations between radar echoes and visual aurorae have not been 
found. Mr Piggott remarked that it was unfortunate that there was not time to 
deal with magnetic and ionospheric disturbances following nuclear explosions in 
the upper atmosphere where the source of disturbance was known. In analysing 
IGY data from different stations it should be remembered that the ionosphere 
behaves differently at each place: account should be taken of the known local 
properties and the relation to the geomagnetic poles. 


A Geophysical Discussion on the Geodetic Uses of Artificial Satellites 
was held on 1961 March 24, under the chairmanship of Dr A. H. Cook. In open- 


ing the meeting the Chairman said that the discussion was concerned with the 
ways in which geometrical relations between points on the surface of the Earth 
could be derived from the observation of satellites. The discussion would not 
consider the information which could be obtained about the external gravity field 
of the Earth from the motions of satellites; a subject which had been discussed 
on a number of occasions in the Society. The fundamental geodetic use of satellites 
is their simultaneous observation from stations whose relative positions we wish 
to determine. In this way one may get information about the shape of the surface 
of the Earth and of its external equipotential surfaces, which could not be obtained 
from ordinary geodetic methods confined as they were to land areas. So far no 
satellites have been used for these geodetic purposes. 

The first speaker, Brigadier G. Bomford, outlined the objectives of the work 
and estimated the accuracy with which the observations must be made. ‘The 
major problem is the connection of the reference spheroid across the oceans. 
Thus for example the form of the American and European arbitrary spheroids is 
well known, but they are not on the same reference datum. Radar links already 
exist across the Atlantic, but these only link latitude and longitude leaving the 
distances to the centre of the spheroids indeterminate. In order to reduce two 
surveyed continents to the same reference datum, simultaneous observations 
are made of a satellite from stations on the two continents and hence the 
geometrical relations between the observing stations obtained. The main 
alternative to this method is the indirect method of observing gravity which until 
very recently has in any case been impracticable in surface ships. We know that 
the geoid departs from any spheroid by + 50m, possibly + 100 m, but not by more. 
Since it is these small scale features which we wish to observe, accuracy to within 
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I m may at present be regarded as perfection and to within 10 m adequate, while an 
observation which might be 30 to 40m in error would not be worth very great 
expense. From the geophysical point of view the slope of the geoid is also of inter- 
est, while cartographically the determination of latitude and longitude is important, 
but for both of these investigations standards can be less stringent. Thus for the 
relative continental links horizontal fixes to within 30m would be valuable. No 
matter what the final aim of the observation, the actual intersection achieved by 
stations involves a series of simultaneous photographs against a stellar background. 
For errors at 2000km with an orthogonal fix an accuracy of 1” is essential while 
0”-2 would be desirable. It appears that the timing must be by the satellite’s own 
flash, since with a velocity of 7000 m/s it cannot otherwise be achieved with the 
necessary precision. After mentioning the interest in the Mediterranean Crete— 
N. Africa link and the far more ambitious projects in the Far East and the Southern 
Pacific the speaker concluded by drawing attention to the need for detailed ground 
administration. 

Dr Atkinson began his contribution by describing how in the simplest optical 
case of a flashing satellite, the flash is photographed against the stars by cameras 
moving equatorially so that both the flash and the stars appear as points. Timing 
is then unnecessary although it can add further information: gravity anomalies 
are totally irrelevant and refraction almost so (Ann.I.G.Y. XII Pt I, 201.). 
Two “known” stations suffice to locate a flash, and two flashes suffice to locate 
any unknown stations that observes them both. The practical upper limit to the 
zenith distance is about 70°, and the Table gives the height h, slant range r, and 
flash-power p, in terms of the ground-distance, D, from the sub-flash point to the 
observer, all for this Z.D. 


D(km) 500 1000 1500 2000 2500 3000 
h(km) 208 474 808 1225 1747 2404 
r(km) 549 1139 1782 2496 3304 4236 
P(c.sec.) 134 576 1411 2769 4818 7997 


(The flash-power, in candle-seconds, is computed for an aperture of 1m and 
ilford S.R. emulsion.) Satisfactory intersection-angles can be obtained with a 
total distance from each known station to the unknown one of about 1-7D; each 
flash is then distant D from the unknown station and from one known one, and is 
considerably nearer the other known one, at a fairly small Z.D. 

The accuracy is limited by the seeing, since the flash-image is not averaged 
over the whole exposure as the star-images are. Information on the seeing at large 
zenith distances is rather scanty and uncertain; a number of Greenwich plates 
taken at zenith distances up to 66° have, however, been examined, and the smallest 
measurable star-images have diameters of 4"-2 (+ 0"-6); this extrapolates to about 
4"°8 at 70°. Estimating the p.e. of the random displacements which produce such 
a disk as about 0-8 of the disk-radius, and allowing a little off the observed radius 
for diffusion in the emulsion, we get + 1”-8 for the p.e., in each co-ordinate, of any 
observation at 70°. 

If we take D = 2000km we have r = 2500; +1°°8 at 2500km is +22m in 
the directions perpendicular to this ray. The ray from the other known station is 
shorter and has better seeing, and the p.e. here is say + 10m. Two flashes, each 
located in this way, will locate an unkown station with p.e. of +30 to + 40m, de- 
pending on the intersection angles; + 40m at 1-7 x 2000km is one part in 85 000. 
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Multiple flashes, and an additional “known’’ station using smaller Z.D.’s only, 
can then give one in 2 or 3 x 10°, and 3 400km is enough for a direct link between 
West Africa and Brazil. 

If observations are entirely confined to smaller Z.D.’s the accuracy improves 
considerably. This will happen if high satellites are used for short lines; it will 
also happen if orbital information can be used on long lines, as may be possible 
if two stations at each end of a long line measure the absolute speed of the satellite 
when near them and the mean speed is then derived from these terminal speeds 
and from orbital data. Non-flashing satellites, if bright enough (Echo I) can be 
used in all these techniques, the trail being repeatedly interrupted at known times, 
e.g. by a rotating shutter with a millisecond chronograph. Long east-west lines 
cannot, however, be obtained, since all stations must be in darkness and the 
satellite in sunlight. 

Dr A. P. Willmore then spoke on the technique of optical observations. Since 
the satellite when viewed from the Earth’s surface moves with an angular velocity 
which is neither constant nor even accurately predictable, the problem is to obtain 
a sufficient exposure to record its position. The sensitivity of the optical detecting 
systems may be limited in any one of the following three ways: by the necessity of 
recording at least a few photon events in each image element, by statistical fluctua- 
tions in the background due to the night sky noise, or by system noise (chemical 
fog, for example, in a photographic emulsion). When these limitations are exam- 
ined for any realistically large, good quality optical system, it is found that at high 
quantum yields characteristic of a photoelectric detector, the threshold sensitivity 
is always determined by the night sky background, whilst at the lower quantum 
yield of the photographic emulsion, chemical fog is the critical factor. The prac- 
tical limits of detection are g™ for photographic and 11™ for photoelectric detector. 
By tracking the camera at the predicted angular velocity of the satellite a further 
improvement of 4~5™ can be obtained. The only precise tracking system of high 
sensitivity in regular operation is photographic. It is operated by the Smith- 
sonian Astrophysical Observatory using specially designed Baker-Nunn cameras. 
These have an f/1 super-Schmidt optical system of 20-in. focal length. ‘The cameras 
are tracked at the predicted angular rate and satellites of 13™ have been recorded. 
However systems employing photoelectric detection, which in addition to their 
higher sensitivity lend themselves readily to very accurate timings, are under 
active study—one for example, using image intensification, by Professor Magee 
at Imperial College, London and another using photomultipliers at University 
College, London. 

Dr Groves gave an account of the principles and the application of the Transit 
satellite. Measurements of the change in the observed frequency of a radio 
transmitter carried in the satellite during transit have already been undertaken for 
most of the satellites which have been launched. Most of the frequency change 
arises from the doppler effect and an analysis of this frequency change enables the 
time and distance of nearest approach to be obtained. Contributions to the fre- 
quency change are however present on account of ionospheric refraction and 
oscillator drift, and these are limitations on the accuracy with which the satellite’s 
orbit can be determined. The ionospheric effect is corrected with information 
obtained by transmitting on two or more frequencies. Above 100 Mc/s 2nd and 
higher order terms in the expression for frequency shift due to refraction may be 
neglected so that by transmitting on two frequencies the correction is estimated 
to be accurate to within o-1 c/s. To improve oscillator stability special attention 
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Im may at present be regarded as perfection and to within 10m adequate, while an 
observation which might be 30 to 40m in error would not be worth very great 
expense. From the geophysical point of view the slope of the geoid is also of inter- 
est, while cartographically the determination of latitude and longitude is important, 
but for botk of these investigations standards can be less stringent. Thus for the 
relative continental links horizontal fixes to within 30m would be valuable. No 
matter what the final aim of the observation, the actual intersection achieved by 
stations involves a series of simultaneous photographs against a stellar background. 
For errors at 2000km with an orthogonal fix an accuracy of 1” is essential while 
o”-2 would be desirable. It appears that the timing must be by the satellite’s own 
flash, since with a velocity of 7000 m/s it cannot otherwise be achieved with the 
necessary precision. After mentioning the interest in the Mediterranean Crete— 
N. Africa link and the far more ambitious projects in the Far East and the Southern 
Pacific the speaker concluded by drawing attention to the need for detailed ground 
administration. 

Dr Atkinson began his contribution by describing how in the simplest optical 
case of a flashing satellite, the flash is photographed against the stars by cameras 
moving equatorially so that both the flash and the stars appear as points. Timing 
is then unnecessary although it can add further information: gravity anomalies 
are totally irrelevant and refraction almost so (Ann.I.G.Y. XII Pt I, 201.). 
Two “known” stations suffice to locate a flash, and two flashes suffice to locate 
any unknown stations that observes them both. The practical upper limit to the 
zenith distance is about 70°, and the Table gives the height h, slant range r, and 
flash-power p, in terms of the ground-distance, D, from the sub-flash point to the 
observer, all for this Z.D. 


D(km) 500 1000 1500 2000 2500 3000 
h(km) 208 474 808 1225 1747 2404 
r(km) 549 1139 1782 2496 3304 4236 
P(c.sec.) 134 576 1411 2769 4818 7997 


(The flash-power, in candle-seconds, is computed for an aperture of 1m and 
Ilford S.R. emulsion.) Satisfactory intersection-angles can be obtained with a 
total distance from each known station to the unknown one of about 1°7D; each 
flash is then distant D from the unknown station and from one known one, and is 
considerably nearer the other known one, at a fairly small Z.D. 

The accuracy is limited by the seeing, since the flash-image is not averaged 
over the whole exposure as the star-images are. Information on the seeing at large 
zenith distances is rather scanty and uncertain; a number of Greenwich plates 
taken at zenith distances up to 66° have, however, been examined, and the smallest 
measurable star-images have diameters of 4”-2 (+ 0"-6); this extrapolates to about 
4° ‘8 at 70°. Estimating the p.e. of the random displacements which produce such 
a disk as about 0°8 of the disk-radius, and allowing a little off the observed radius 
for diffusion in the emulsion, we get + 1"-8 for the p.e., in each co-ordinate, of any 
observation at 70°. 

If we take D = 2000km we have r = 2500; +1”°8 at 2500km is +22m in 
the directions perpendicular to this ray. The ray from the other known station is 
shorter and has better seeing, and the p.e. here is say + 10m. Two flashes, each 
located in this way, will locate an unkown station with p.e. of +30 to + 40m, de- 
pending on the intersection angles; + 40m at 1-7 x 2000km is one part in 85 000. 
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Multiple flashes, and an additional “known” station using smaller Z.D.’s only, 
can then give one in 2 or 3 x 10°, and 3 400km is enough for a direct link between 
West Africa and Brazil. 

If observations are entirely confined to smaller Z.D.’s the accuracy improves 
considerably. This will happen if high satellites are used for short lines; it will 
also happen if orbital information can be used on long lines, as may be possible 
if two stations at each end of a long line measure the absolute speed of the satellite 
when near them and the mean speed is then derived from these terminal speeds 
and from orbital data. Non-flashing satellites, if bright enough (Echo I) can be 
used in all these techniques, the trail being repeatedly interrupted at known times, 
e.g. by a rotating shutter with a millisecond chronograph. Long east-west lines 
cannot, however, be obtained, since all stations must be in darkness and the 
satellite in sunlight. 

Dr A. P. Willmore then spoke on the technique of optical observations. Since 
the satellite when viewed from the Earth’s surface moves with an angular velocity 
which is neither constant nor even accurately predictable, the problem is to obtain 
a sufficient exposure to record its position. The sensitivity of the optical detecting 
systems may be limited in any one of the following three ways: by the necessity of 
recording at least a few photon events in each image element, by statistical fluctua- 
tions in the background due to the night sky noise, or by system noise (chemical 
fog, for example, in a photographic emulsion). When these limitations are exam- 
ined for any realistically large, good quality optical system, it is found that at high 
quantum yields characteristic of a photoelectric detector, the threshold sensitivity 
is always determined by the night sky background, whilst at the lower quantum 
yield of the photographic emulsion, chemical fog is the critical factor. The prac- 
tical limits of detection are g™ for photographic and 11™ for photoelectric detector. 
By tracking the camera at the predicted angular velocity of the satellite a further 
improvement of 4-5™ can be obtained. The only precise tracking system of high 
sensitivity in regular operation is photographic. It is operated by the Smith- 
sonian Astrophysical Observatory using specially designed Baker-Nunn cameras. 
These have an f/1 super-Schmidt optical system of 20-in. focal length. The cameras 
are tracked at the predicted angular rate and satellites of 13™ have been recorded. 
However systems employing photoelectric detection, which in addition to their 
higher sensitivity lend themselves readily to very accurate timings, are under 
active study—one for example, using image intensification, by Professor Magee 
at Imperial College, London and another using photomultipliers at University 
College, London. 

Dr Groves gave an account of the principles and the application of the Transit 
satellite. Measurements of the change in the observed frequency of a radio 
transmitter carried in the satellite during transit have already been undertaken for 
most of the satellites which have been launched. Most of the frequency change 
arises from the doppler effect and an analysis of this frequency change enables the 
time and distance of nearest approach to be obtained. Contributions to the fre- 
quency change are however present on account of ionospheric refraction and 
oscillator drift, and these are limitations on the accuracy with which the satellite’s 
orbit can be determined. The ionospheric effect is corrected with information 
obtained by transmitting on two or more frequencies. Above 100 Mc/s 2nd and 
higher order terms in the expression for frequency shift due to refraction may be 
neglected so that by transmitting on two frequencies the correction is estimated 
to be accurate to within o-1 c/s. To improve oscillator stability special attention 
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needs to be given to thermal control. A stability of 1 in 10°, or 0-2c/s at 200 Mc/s 
is not far short of being realized by rocket-borne instrumentation. The Transit 
satellites have been designed to overcome the difficulties of oscillator drift and 
ionospheric refraction so that an accurate orbit may be computed from a network 
of tracking stations (Guiep, W. H., Weiffenvach, G. C., Proc. Inst. Rad. Eng. 48, 
507, 1960). The orbit may then in turn be used for locating the position of an 
observer. Observations from six doppler stations were used to derive the orbit 
of Transit Ib by combining a single day’s observations to determine the initial 
position and velocity of the satellite for that day, the orbit being derived by numeri- 
cal integration (Cohen C. J. & Anderle R. J., Science, 132, 807, 1960). It was 
found necessary to take account of the third harmonic of the Earth’s gravitational 
potential and the corresponding accuracy of the orbit appears to be of the order 
of }km. In using the Transit satellites for the location of the position of observers 
it is important that the satellite give sufficient coverage, have a useful life and 
avoid local anomalous drag effects. For all these reasons a height of at least 800 km 
is desirable. Ten stations would be required to give the data for the computation 
of the orbit which could be transmitted twice a day and hence a record of predicted 
positions of the satellite given every twelve hours. To use the satellite for navi- 
gation the observer would need to be able to digitize the doppler information and 
use a computer to obtain his position knowing the orbit and using his doppler 
information. The speaker finally investigated the effect of systematic errors in 
frequency on position determination. Assuming the orbit to be accurate, it is 
found that a doppler record extending over several minutes would enable the 
observer's position to be determined to approximately 1oom if the frequency 
variation is 1 in 10%. The accuracy with which position along a line of latitude is 
obtained would be degraded by a factor of 10 in the case of a polar orbit passing 
overhead. 

The Chairman regretted that there was so little time for general discussion and 
asked Dr R. L. Waterfield to give his contribution on the problem of obtaining 
accurate timed positions of satellites. Using the technique of interrupting the 
photographic trail and timing the interruptions an accuracy of about 7” has been 
achieved, most of which was due to vibrations set up by the heavy shutter. Subse- 
quently this vibration iias been considerably reduced so that accuracies of 5” are 
possible. By making the moving parts of the shutter of much lighter material the 
error should approach 1”. With the direct recording of the Rugby time signal on 
the chronograph it is hoped to measure the “‘photographic centre of gravity” of 
the short measured sequences to better than 0-01 s. Dr Waterfield thought that 
the large number of measured segments and the length of each one (1’ or more) 
should make thé error due to bad seeing, referred to by Dr Atkinson, negligible 
even at low altitudes. Using a 6” Cooke triplet of 26” focus, with an Oa Kodak 
emulsion an object with the angular velocity of Echo should reach about magni- 
tude 6, while the faster moving objects might reach magnitude 4. 

Mr B. C. Browne was then asked to comment on the requirements from satellites 
of geophysicists working at sea. He said that for geophysical work at sea, an 
accurate knowledge of position is important. If gravity measurements are being 
made, it is also necessary to know one’s east-west velocity, relative to the Earth, to 
within a fraction of a knot. The “Transit” system, which promises to give a fix 
within half a mile every go minutes, is therefore most attractive. Preliminary 
results for Transit I are encouraging, but some observations made at Cambridge 
suggest that the predicted time of closest approach has a small systematic error. 
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Could this be due to some uncertainty in the longitude of the United States 
relative to Europe? 

Closing the discussion, the Chairman said that it seemed to him that three 
points needed a great deal more attention before it would be possible to make 
geodetic observations by optical methods. A detailed study should be made of 
the orbits that could be used, bearing in mind the visibility of the satellite and the 
places between which connexions were required; the requirements of the ground 
organization should be formulated in detail; and the problems of a world-wide time 
service giving times correct to 1 ms should be examined. 
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Book Review 


Report on the 
Symposium on July 1959 Events and Associated Phenomena 


(Monograph No. 7, 1960 November, published by the I1.U.G.G.) 


The Symposium on the 1959 July Events and Associated Phenomena was held 
as part of the scientific programme of the International Association of Geo- 
magnetism and Aeronomy during the 12th General Assembly of the I.U.G.G. 
The suggestion for an interdisciplinary study of the relationships between various 
scientific observations in the field of solar—terrestrial relationships was made by 
C.O.S.P.A.R., which called attention to the need for a detailed study of some 
significant solar—terrestrial events. Accordingly, Working Group 2 of C.O.S.P.A.R. 
decided to select the month of July 1959 for a special study of this kind. The 
reason for this selection was not only the occurrence of three magnetic storms 
accompanied by cosmic ray intensity decreases, but also the occurrence of polar 
cap absorption preceding these three storms and the observation of additional low 
energy cosmic radiation at balloon altitudes. 

The sequence of events began with three solar flares of importance 3+ com- 
mencing some time shortly before 0210 hours U.T. on July roth, at 0319 hours 
U.T. on July 14th, and at 2118 hours U.T. on July 16th. These flares were 
followed by three distinct geomagnetic storms with sudden commencements at 
1625 hours on July 11th, at 0803 hours on July 15th, and at 1638 hours on July 
17th. There was some discussion at the symposium on the possibility that a 
sharp change in the length of the day occurred at this time, but the evidence for 
this seems to be rather doubtful. 

Each of the three magnetic storms was accompanied by a cosmic ray decrease 
of the Forbush type, the third of these resulting in the lowest cosmic ray intensity 
recorded during the present solar cycle. Solar protons were observed by balloon- 
borne cosmic ray recorders in association with each of the three flares, but follow- 
ing the flare of the 17th a small increase in cosmic ray intensity was recorded by 
high latitude stations at sea level indicating the presence on this occasion of solar 
protons with momenta up to 1Gev/c. Carmichael & Steljes at Chalk River, 
Ontario, have assembled a very complete collection of cosmic ray data from many 
stations all over the world and these are available in the form of photographic 
transparencies. 

The Symposium included a number of papers on the fine structure of the mag- 
netic storms, the aurora and the changes in ionospheric characteristics during this 
period, together with an account of solar emission at radio frequencies and a 
summary of the solar data generally for this period. This published account of 
the meeting contains a store of related data covering this period of remarkable 
solar activity which will be extremely useful to all who are interested in the 
problems of solar—terrestrial relations. 

H. ELLioT 
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ERRATUM 


Geophys. J., 3, No. 4, 1960: 
Report of the XII General Assembly of the International Union of Geodesy & 
Geophysics. Section 4, Age Determination, p. 473, lines 26 and 30: 


For ‘““Thorium-—231 (Protoactinium)”’ read 
“Protoactinium-231” 
and for “Th®!” read 
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Barrett and C. Batey (Oxford University Press, 1954). In particular the abbrevia- 
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tracing film. Original drawings must not be lettered; lettering should be indicated 
on copies or photoprints which must be supplied in addition to the originals. 
Wherever possible, lettering should be kept outside the diegram so that it may b: 
set in type. Diagrams should be drawn at twice to three times the size at which 
they will be printed. The maximum dimensions for diagrams, when reduced and 
with their lettering and captions, are 8 in by 5 in. A wide variety of flat tints of 
dots, lines and shadings can be applied to line drawings by the blockmaker, Photo- 
graphs for reproduction should be unmounted giossy prints and should be accom- 
panied by lettered prints. 


5. References should be quoted in the form recommended by the Royal Society 
but authors may also quote titles of papers if they wish. References in the text are 
made in the form (Smith 1940). 


6. Authors whose work involves much mathematics should follow carefully the 
recommendations in The Printing of Mathematics. 

7. The Editors will be glad to advise authors on any special points of difficulty 
arising in the prepcration of papers for the Journal. 
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